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ABSTRACT 


This  investigation  is  concerned  with  miniTnnm  weight  designs  of  arch 
structures.  Using  the  finite  element  method,  the  arch  is  modeled  by 
contiguous  bar-beam  elements.  Element  stiffness  coefficients  in  terms  of 
local  degrees  of  freedom  are  related  to  system  degrees  of  freedom  through 
local  to  global  coordinate  transformations.  After  coordinate  transformations, 
element  stifffiess  coefficients  are  assembled  into  FEM  stiffiiess  equations  for 
the  arch  structtire.  An  objective  function  for  weight  minimization,  with 
constraints  on  failure,  arch  geometry,  and  section  dimensions,  is  minimized 
by  the  DOT  optimization  code.  Results  are  presented  for  a  number  of  cases. 
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1.  INTRODUCTION 


The  arch  has  been  employed  in  the  fabrication  of  engineering 
structures  for  over  five  thousand  years.  Its  suitability  to  compressive  load 
design  made  it  especially  favored  when  masonry,  not  steel,  was  the  principle 
building  material.  (Due  to  masoniy  being  stronger  in  compression  than  in 
tension.)  Its  elegant  shape,  more  natural  and  graceful  than  the  straight  lines 
and  perpendiculars  of  traditional  man-made  structures,  made  it  fashionable 
among  architects  and  civil  engineers.  Its  influence  can  be  dted  in  diverse 
cultures,  among  which  include  the  Egyptians,  Mesopotamians,  Romans,  (see 
Figure  1.1),  Byzantines,  French,  Chinese,  and  English. 

A  number  of  investigations  on  the  optimization  of  arches  have  been 
conducted  over  the  years.  In  1976,  Farshad  [Ref.  1],  using  the  calculus  of 
variations,  derived  optimality  conditions  in  the  form  of  nonlinear  partial  dif¬ 
ferential  equations  for  hinged-hinged  arches.  An  augmented  functional, 
comprised  of  the  total  potential  energy  of  the  system  and  the  objective 
function,  appended  to  the  system  via  Lagrange  multipliers,  when  minimized 
with  respect  to  state  variables  and  with  respect  to  design  variables  yield  the 
system  equilibrium  equations,  and  the  optimality  conditions  respectively. 
Three  objective  functions  were  imposed: 

-  optimal  thrust 

-  minimum  length  of  arch 
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minimum  volume 


The  arch  span  and  the  loading  are  specified.  The  nonlinear  system  of 
optimality  equations  were  presented  but  not  solved. 

In  1980,  Rozvany  et  al  [Ref.  2]  considered  the  problem  of  arch 
optimization  using  the  Prager-Shield  criteria.  Here,  the  arch  was  in  fact  a 
funicular  frame  with  beams  rigidly  interconnected  to  one  another.  Only 
statically  determinate  systems  were  investigated.  The  first  "arch"  with  a 
specified  span  consisted  of  two  inclined  beams  with  a  concentrated  load  along 
the  center  of  symmetry.  The  second  investigation  dealt  with  an  "arch" 
consisting  of  three  beam  segments,  the  center  segment  being  horizontal,  and 
inclined  members  from  the  hinged  supports.  Two  concentrated  loads  were 
applied  at  the  intersections  of  the  inclined  members  with  the  horizontal 
center  member.  For  the  single  load  "arch"  it  was  foimd  that  the  optimal 
"arch"  develops  either  bending  only  or  axial  forces  only  in  the  entire  structure 
depending  on  the  range  of  the  4Lyd  ratio,  where  L  is  the  span  of  the  structure 
and  d  is  the  constant  depth  of  the  cross-section.  For  4Li/d  greater  than  8,  the 
optimal  structure  has  a  height  half  of  the  span,  and  there  is  only  axial  force 
throughout  the  structure.  For  4iyd  less  than  8,  the  optimal  structure  is  a 
straight  horizontal  beam  (i.e.,  the  height  is  zero),  and  only  bending 
throughout  the  structure.  The  width  of  the  beam  segments  for  the  optimal 
"arch"  varies  linearly  fimm  the  hinged  support  to  the  center  line. 
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In  a  1980  paper,  Lipson  et  al  [Ref.  3]  investigated  the  optimal  design  of 
arches  using  the  complex  method.  Both  the  arch  shape  and  the  cross 
sectional  dimensions  are  the  design  variables  for  the  minimum  weight 
structure.  Only  symmetric  arches  with  constant  depth  and  constant  width 
were  considered.  The  section  was  taken  as  a  thin  walled  rectangular  tube 
with  vertical  and  horizontal  wall  thicknesses  as  design  variables.  The  arch 
was  approximated  with  equal  length  straight  beam  sections.  Thus,  each 
beam  segment  had  three  design  variables,  the  two  wall  thicknesses  and  the 
left  end  vertical  location.  In  addition,  the  uniform  height  and  uniform  width 
of  the  rectangular  tube  were  design  variables.  Side  constraints  in  the  form  of 
upper  and  lower  bounds  were  placed  on  all  of  the  design  variables.  A 
modified  version  of  the  complex  method  of  Box  was  used  as  the  scheme  to 
obtain  a  "fully-stressed"  optimum  design.  The  shape  of  the  arch  was  taken  as 
a  parabola.  The  optimization  algorithm  provided  the  minimum  weight 
parabolic  arch  for  a  uniform  load  over  a  specified  span.  It  was  shown  that  a 
parabolic  arch  with  a  rise  0.342  times  the  span  length  is  the  optimal 
parabolic  shape  for  the  case  of  a  uniform  horizontal  load.  Deviations  within 
10%  from  this  rise  have  negligible  effect  on  the  optimal  weight.  It  was 
further  shown  that  parabolic  steel  arches  wiU  fail  due  to  their  own  weight  at 
span  lengths  greater  than  1,543  ft.  For  relatively  high  arches,  the  maximum 
axial  thrust,  which  occurs  at  the  supports,  approaches  half  of  the  total 
uniform  load.  The  results  of  a  parametric  study  of  optimal  steel  arches  are 
presented. 

In  a  1988  paper,  Ang  et  al  [Ref.  4]  investigated  the  optimal  shape  of  an 
arch  under  bending  and  axial  compression.  The  cross-section  of  the  arch  was 
rectangiilar  with  specified  constant  depth  and  variable  width.  With  the 
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centroidal  axis  of  the  arch  given  by  the  summation  of  products  of  cubic  spline 
functions  vdth  linear  coefficients,  the  arch  axis  is  permitted  to  take  on  any 
smooth  shape.  The  linear  coefficients  in  these  products  are  design  variables 
to  be  determined  in  the  optimization  process.  The  authors  considered  arches 
imder  a  uniformly  distributed  horizontal  load  with  three  types  of  boundaries, 

-  simply  supported-simply  supported 

-  clamped-damped 

-  clamped-simply  supported 

A  yield  failure  constraint  was  imposed.  A  new  technique  for  smoothing  the 
objective  function  is  presented.  It  was  shown  that  the  optimal  shape  of  the 
arch  is  a  parabola  with  a  rise  equal  to  0.433  time  the  span  of  the  arch.  No 
other  results  are  presented.  It  should  be  noted  that  the  results  of  [Ref.  3]  and 
[Ref.  4],  with  regard  to  the  ratio  of  rise  to  span  for  an  optimal  parabolic  arch, 
do  not  agree. 

In  the  study  of  arches,  it  is  necessary  to  determine  how  they  will  be 
defined.  One  prevalent  school  of  thought  defines  an  arch  as  a  curved 
structure,  which  when  supported  at  both  ends  and  loaded  vertically  develops 
horizontal  reactions.  This  is  apparently  intended  to  eliminate  thick  walled 
curved  beams  and  straight  beams  which  develop  (virtually)  no  horizontal 
reactions  when  loaded  laterally. 

Others  define  an  arch  as  a  curved  beam  whose  cross-sectional 
dimensions  are  small  relative  to  its  radius  of  curvature.  Hence,  the 
centroidal  and  neutral  axes  are  assumed  to  coincide.  How  the  structure  is 
loaded  and  supported  becomes  secondary.  This  description  was  chosen  to 
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facilitute  the  development  of  a  finite  element  code  capable  of  generating 
horizontal  and  vertical  displacements  and  slopes  for  an  arbitrarily  loaded 
arch.  Without  the  thin  depth  assumption,  complications  arise  in  the 
calculations  of  the  slopes  and  displacements  because  the  arch  will  no  longer 
behave  as  predicted  by  the  beam  equilibrium  equation: 

(Elv"r=py(s)  (1.1) 

and  the  bar  equilibrium  equation: 

(AEuO  =  ~-p,(s)  (1.2) 

where  E  represents  Yotmg^s  Modulus,  I  the  cross-sectional  moment  of  inertia, 
A  the  cross-sectional  area,  v  the  lateral  displacement,  u  the  axial 
displacement,  p^  the  axial  loading,  and  p^  the  lateral  loading. 

Once  the  displacements  and  slopes  are  determined,  the  local  stresses 
can  be  calculated  throughout  the  member  using  appropriate  stress- 
displacement  relations.  Thus  able  to  determine  the  stress  distribution,  the 
arch  may  be  designed  to  achieve  minimum  volume  (and  hence  weight)  while 
maintaining  the  developed  stresses  below  some  predefined  maximum 
allowable  stress  value. 

The  aim  of  this  study  is  to  "optimize"  a  linear,  elastic,  isotropic,  and 
homogeneous  arch  imder  a  variety  of  loadings  and  end  conditions.  Although 
these  limitations  are  not  physically  essential,  they  were  necessaiy  to  make 
the  investigation  tenable  given  the  time  constraints  of  thesis  research 
activity.  Optimization  in  this  investigation  will  refer  to  the  variance  of  the 
cross-sectional  geometry  to  achieve  a  more  uniform  stress  distribution 
throughout  the  member.  This  results  in  less  material  used  and  a  more 
efficient  structure.  VMA  Engineering's  Design  Optimization  Tools  [Ref.  5]  is 
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used  to  perform  the  optimization  subject  to  prescribed  constraints  on  the 
design  variables  as  well  as  stress  limitations.  The  objective  function  to  be 
minimized  is  the  total  volume  of  the  arch  while  maintaining  stresses  below 
the  yield  strength  of  the  arch  material.  Though  a  rather  simplistic  model,  it 
forms  a  foundation  upon  which  further  research  into  more  complex 
geometries  and  conditions  may  be  developed. 
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11.  PROBLEM  FORMULATION 


Perhaps  the  most  common  optimization  in  structural  mechanics  is  the 
minimization  of  an  element's  weight,  subject  to  a  specified  loading.  Such  will 
be  the  case  for  this  investigation.  To  make  the  investigation  tenable,  the 
problem  needed  to  be  narrowed  down  in  its  scope.  The  assumptions  and 
approximations  made  in  this  study  are: 

-  The  arch  is  approximated  by  a  series  of  straight  bar/beam  elements 
which  behave  according  to  the  beam  equation  (1.1)  and  bar  equation 
(1.2).  (See  Figure  2.1) 

-  The  arch  material  is  isotropic,  homogeneous,  and  linearly  elastic. 

-  The  arch's  cross-sectional  area  will  always  be  of  a  solid  rectangular 
geometry. 

-  The  arch  has  a  constant  circular  radius  of  curvature. 

-  The  arch  "fails"  if  its  internal  stresses  exceed  the  yield  strength,  S^. 

It  should  be  noted  that  the  third  and  fourth  assumptions  are  not  inherent  to 
the  general  optimization  problem  but  rather  are  imposed  to  limit  the  scope  of 
this  initial  investigation.  Follow  on  investigations  will  relax  these 
restrictions. 
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Figure  2.1:  Bar-Beem  Model  of  the  Arch 


Althoiigh  these  suppositions  limit  the  applications  for  which  the  optimization 
can  be  used,  they  form  a  foundation  upon  which  further  research  can  be 
based. 


A.  THE  OBJECTIVE  FUNCTION 

The  objective  of  this  investigation  is  the  minimization  of  the  structure's 
weight  while  maintaining  a  stress  distribution  which  does  not  exceed  the 
yield  strength.  Since  the  weight  of  the  arch  is  directly  proportional  to  the 
volume  of  material  from  which  it  is  made,  the  objective  will  be  satisfied  if  the 
total  volume  of  the  arch  is  minimized.  That  is: 

Objective  *  MmSbihil^  (2.1) 

where  b^,  h^,  and  1^  represent  the  width,  height,  and  length  of  the  i*^  element, 
respectively,  and  NEL  is  the  total  number  of  elements.  (See  Figure  2.1b) 
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With  an  objective  function  defined,  the  next  step  in  the  optimization 
process  is  to  impose  any  necessary  (design)  constraints  upon  the  system.  The 
constraints  m\ist  be  provided  to  represent  the  assumptions  contained  in  the 
modeling  equations.  They  should  be  utilized  to  avoid  undesirable  behavior 
such  as  buckling  and  yielding.  They  may  also  be  used  to  apply  any 
limitations  on  behavior  as  desired  by  the  designer.  The  constraints  employed 
in  this  investigation  follow. 

B.  THE  STRENGTH  CONSTRAINT 

This  study  assumes  a  linearly  elastic  arch.  Therefore,  the  applied 
loading  must  not  cause  the  structure  to  exceed  the  elastic  limit  of  the 
material  from  which  it  is  made.  Hence,  a  design  constraint  which  will 
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prevent  this  mode  of  "failure  by  yielding"  must  be  imposed.  To  fulfill  this 
constraint,  the  internal  stresses  developed  must  remain  below  the  3rield 
strength  of  the  material.  Defining  the  maximum  stress  developed  in  the  i‘** 
element  of  the  arch  to  be  and  the  yield  strength  of  the  arch  material  as  S^, 

this  constraint  becomes 

a-  <  S 

V/j  — 

or  in  dimensionless  form: 

o./Sy-1^0  (2.2) 

C.  GEOMETRIC  CONSTRAINTS 

To  use  the  beam  and  bar  equations,  limits  must  be  imposed  upon  the 
geometric  dimensions  of  the  structure.  It  is  compulsory  that  the  depth  and 
width  be  of  at  least  an  order  of  magnitude  smaller  than  the  radius  of 
curvature  for  the  beam  and  bar  equations  to  be  applicable.  Hence, 

hj  ^  Rj/lO 

or  in  dimensionless  form: 

lOhj/Rj-l^O  (2.3) 

Similarly,  as  the  width  increases,  the  arch  approaches  the  geometry  of 
a  shell.  To  maintain  the  geometry  of  an  arch,  an  imposition  upon  the  width 
dimension  is  also  necessary.  To  avoid  shell  behavior,  a  third  constraint  is 
imposed 

b,<3h, 

or  in  dimensionless  form: 

b./3hj-1^0  (2.4) 
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D.  SIDE  CONSTRADJTS 


Finally,  side  constraints  must  be  assigned  to  the  dimensions  of  the 
design  variables.  For  the  sake  of  simphdty,  this  investigation  will  only  take 
up  the  variation  of  the  cross-sectional  width  dimension  b-.  The  side 

constraints  must  ensure  real  solutions  are  obtained,  i.e.,  the  arch  is  a 
physical  object  and  therefore  h^  and  bj  cannot  be  less  than  a  realistic  finite 

value.  The  side  constraints  chosen  to  reflect  these  limitations  include: 


0.03  in.  $  h^ 

(2.5) 

0.03  in.  ^  b{ 

(2.6) 

Other  constraints  could  have  also  been  considered.  Global  buckling 
and  local  crippling  are  to  name  but  two.  However,  the  cases  to  be  studied  do 
not  warrant  such  a  thorough  delineation.  Therefore,  the  design  and  side  con¬ 
straints  have  been  limited  to  those  dted. 

E.  OPTIMIZAnON  SOFTWARE 

With  a  multitude  of  preprogrammed  optimization  routines  available, 
the  Design  Optimization  Tbols  (DOT)  software  was  chosen.  Its  selection  was 
based  upon  availability,  ease  of  use,  and  reputation.  DOT  is  a  FORTRAN  77 
optimization  software  package  available  from  VMA  Engineering.  To  perform 
a  variety  of  optimization  tasks,  DOT  uses: 

-  The  Modified  Method  of  Feasible  Directions, 

-  Broydon-Fletcher-Goldfarb-Shanno  (BFGS)  Variable  Metric  Method, 

-  Polynomial  Interpolation  with  botmds,  and 

-  Sequential  Linear  Programming  (SLP) 
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A  user  provided  "main"  program  is  used  to  input  the  variables  required  by 
DOT.  DOT  in  turn  calls  a  user  provided  subroutine  which  defines  the 
objective  function,  design  constraints,  and  design  variable  side  constraints. 
DOT  iteratively  evaluates  the  objective  function,  refining  the  design  variables 
imtil  the  optimal  solution  is  obtained. 

The  parameters  used  to  calculate  the  objective  function  and  constraints 
must  be  known  before  any  optimization  can  occur.  The  vanables  from 
equations  (2.1)  through  (2.6)  include: 

-  The  number  of  elements  used,  NEL. 

-  The  arch  radius  of  curvature,  R. 

-  The  height  of  the  i^**  element,  hj. 

-  The  length  of  the  i*^  element,  1^. 

-  The  width  of  the  i*^  element,  bj. 

-  The  yield  strength  of  the  material  firom  which  the  arch  is  made,  Sy 

~  The  stress  at  the  i*^  node,  Oj. 

Of  these  terms,  the  number  of  elements  is  the  choice  of  the  designer.  The 
arch  radius  of  curvature  and  height  are  constant  throughout  the  span  of  the 
arch  and  are  defined  by  the  problem.  For  simplicity,  the  length  of  each 
element  is  made  uniform  such  that: 

\  =  eR/(NEL)  (2.7) 

where  0  represents  the  subtended  arc  of  the  arch.  The  yield  strength  is 
determined  by  the  material  used  to  build  the  arch  and  the  width  is  the  design 
variable  to  be  determined  by  DOT. 
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The  stress  distribution  is  not  so  readily  available.  However,  using  the 
beam  and  bar  equations,  a  finite  element  scheme  can  be  developed  to 
determine  the  arch's  displacements  and  slopes  due  to  a  given  loading. 
Knowing  how  the  displacements  and  slopes  change  throughout  the  arch,  the 
stress  at  a  given  point  can  be  calculated  and  the  optimization  performed. 
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ni.  STRESS  ANAIYSIS 


The  objective  of  this  optimization  is  to  minimize  the  total  weight 
(volume)  of  a  load  bearing  arch  subject  to  specified  constraints.  The  strength 
constraint  requires  that  the  stress  at  any  point  does  not  exceed  the  3rield 
strength  of  the  arch  material.  To  avoid  violating  this  constraint,  the  value  of 
the  stresses  at  any  point  must  be  known.  With  this  requirement,  the  stress 
development  is  pursued. 

The  normal  stress  at  any  point  in  the  arch  is  composed  of  normal 
stresses  due  to  bending  (bending  stress)  and  normal  stresses  due  to  the  axial 
forces  (axial  stress)  acting  upon  the  individual  elements.  Figure  (3.1)  depicts 
these  force  interactions.  The  total  normal  stress  is  the  algebraic  sum  of  these 
components. 

(3.1) 

The  arch  can  also  develop  shear  stresses  due  to  shear  forces.  Due  to  the 
geometric  constraint  defined  by  equation  (2.3),  the  shear  stresses  turn  out  to 
be  insignificant  when  compared  to  the  normal  stresses.  Consequently,  the 
shear  stresses  are  ignored.^ 

To  calculate  the  normal  stress  components,  we  must  first  determine 
how  the  elements  behave.  For  a  simple  straight  beam  element,  the  maximum 
normal  stress  due  to  bending  occurs  at  the  outer  fibers  and  is  given  by 

Oj,  =  Mc/I 


1.  See  Appendix  A  for  a  justification  of  the  omission  of  the  shear  stresses. 
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or  in  terms  of  the  beam  equation: 

Oj,  =  (EIv")c^ 

which  reduces  to: 

Oj,  =  Ecv"  (3.2) 

where  c  is  the  distance  from  the  neutral  axis  to  the  outer  fiber  of  the  beam, 
that  is  csh/2.  See  Figure  (3.1b). 


F 


Ficvra  S.1:  Normal  Streaaea  doe  to  Bending  and  Axial  Force 

In  the  same  manner,  the  normal  stresses  due  to  axial  behavior  can  be 
determined.  The  uniform  normal  stress  due  to  axial  forces,  F,  acting  upon  a 
bar  can  be  defined  by: 

o,*F/A 

or  in  terms  of  the  bar  equation, 

o.*(AEu'yA 


which  reduces  to: 


=  Eu' 


(3.3) 


See  Figure  (3.1c), 

Substituting  equations  (3.2)  and  (3,3)  into  equation  (3.1)  yields  the 
linear  equation 

=  Ecv"+  Eu' 

or  simply 

=  E(cv"  +  uO  (3.4) 

where  v"  and  u'  are  to  be  determined  from  the  solution  of  equations  (1.1)  and 

(1.2). 

From  this  development,  we  see  the  normal  stress  is  a  function  of 
Young's  Modulus,  the  height  of  the  beam,  the  first  derivative  of  the  axial 
displacement,  and  the  second  derivative  of  the  lateral  displacement.  Using 
the  Galerkin  finite  element  method,  approximate  values  of  u'  and  v"  can  be 
determined.  With  these  values,  the  stress  distribution  can  be  calculated 
using  Equation  (3,4)  and  the  optimization  may  then  be  conducted. 
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IV.  FINITE  ELEMENT  ANAIYSIS 


In  order  to  determine  the  stresses  developed  for  a  given  loading,  the 
values  of  u'  and  v"  must  be  determined.  These  derivatives  can  be  found  by 
solving  the  beam  and  bar  equations  using  the  Galerkin  finite  element  method 
(FEM).  The  Galerkin  FEM  is  capable  of  directly  solving  systems  of  linear 
differential  equations  while  preserving  their  synunetry. 

A.  THE  BEAM  EQUATION  DEVELOPMENT 

The  beam  equation  (1.1)  is  a  fourth  order  linear  differential  equation 
requiring  continuity.  Therefore,  a  family  of  cubic  shape  functions  are 
necessary  to  maintain  function  and  slope  continuity.  With  this  in  mind,  the 
Galerkin  FEM  is  performed  on  the  beam  equation.  A  finite  element  method 
is  an  approximation  method  which  transforms  the  differential  equation  of  a 
continuous  system  into  a  system  of  linear  algebraic  equations.  The  method 
begins  by  a  discretization,  that  is  a  partition,  of  the  continuous  domain  into  a 
segmented  domain  of  NEL  elements.  Thereafter,  a  three  step  process  takes 
place. 

The  first  step  is  to  form  an  approximate  solution  V, 

v»v  =  Q'*'y  (4.1) 
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where  v  is  the  exact  solution  in  continuous  space  of  the  beam  equation,  v  is 
the  approximate  solution  in  discrete  space,  is  the  transpose  of  the  column 

vector  of  cubic  shape  functions  which  have  the  Kronecker  delta  property,  and 
Y  is  the  vector  of  coefficients  of  lateral  displacements  and  slopes. 

The  second  step  is  to  form  the  residual  of  the  approximation  where: 

R  =  £(v)-p  (4.2) 

where  £  denotes  the  system  operator,  in  this  case  being  the  beam  equation 
such  that 

£(v)  =  [EI(vy]" 

Substituting  the  beam  equation  (1.1)  and  the  equation  (4.1)  into  equation 

(4.2)  yields: 

/2  =  [EI(QV]"'Py(s)  (4.3) 

The  third  step  is  to  form  the  Galerkin  Equations, 

lQ(fl)ds  =  0  (4.4) 

where  0  represents  a  vector  whose  values  are  zero.  Substituting  equation 

(4.3)  into  equation  (4.4)  yields: 

/DQ[EI(Q'^y)"]”ds  -  JdQp/s)  ds  =  0  (4.5) 

Performing  two  successive  integrations  by  parts  upon  equation  (4.5)  yields: 

eCElCgTy)"]'  Ib  -  Q'EI(Q’'y)-  /DQ"EI(QTy)-ds  -  IdQp,(s)  ds  =  0  (4.6) 

where  B  denotes  the  boundsiry  values  of  the  structure  at  each  end  point. 

Since  the  coefficients  of  the  solution  vector  are  constants,  equation  (4.6) 
may  be  rewritten  as: 

Q[EI(Qf )"]'v  Ig  _  Q-EKeTyv  Ig  +  JnQ"EI(Q’')-d8  v  -  JnQP,(s)  ds  =  0  (4.7) 
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Using  shear  given  by  V=EIv"'  and  moment  given  by  M=EIv",  we  define  the 
boimdary  term  load  vectors  Y  and  M  as 

Y  =  Q[eI(QT)"]v|b  (4.8a) 

and 

M  =  Q'EI(QT)"v  Ig  (4.8b) 

In  addition,  we  define  the  S3rstem  stifibiess  matrix  as 

8®  =  JDQ"EI(Q’')'ds  (4,8c) 

and  the  system  force  vector  as, 

5*  =  /oQPy^®)  ‘is  (4.8d) 

and  upon  substituting  equations  (4.8)  into  equation  (4.7)  we  obtain  the 
system  of  linear  algebraic  equations; 

Y-M  +  ?®y-F*sO  (4.9) 

Moving  the  applied  internal  excitation  and  boundary  terms  to  the  right-hand 
side,  such  that 

g»y  =  F*  +  M-Y  (4.10) 

and  defining  the  load  vector  of  internal  and  external  applied  loads  as 

FB  =  F‘’  +  M-Y  (4.11) 

equation  (4.10)  simplifies  to  the  linear  system; 

5»y  =  FB  (4.12) 

where  the  system  bending  stiffiiess  matrix,  jg  constructed  from  the  union 
of  the  i=l,...4^L  elemental  bending  stifhiess  matrices,  and  the  system 
bending  force  vector  F**  is  formed  from  the  union  of  the  i=l,...,NEL  elemental 
bending  force  vectors  f 
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B.  THE  BAR  EQUATION  DEVELOPMENT 


The  FEM  application  to  the  bar  equation  is  similar  to  that  previously 
conducted  with  the  beam  equation.  The  bar  equation,  however,  is  a  second 
order  linear  difTerential  equation  requiring  C®  continuity.  Hence,  only  a 
family  of  linear  shape  functions  are  necessary  to  maintain  hmction  continuity 
in  the  FEM  development. 

Again,  the  basic  steps  of  the  Galerkin  Method  are  conducted. 

First,  the  approximate  solution  a  is  formed; 

u»a  =  GTu  (4.13) 

where  u  is  the  exact  solution  of  the  bar  equation,  u  is  the  approximate 
solution,  G*^  is  the  transposed  column  of  linear  shape  functions  with  the 

Kronecker  delta  property,  and  u  is  the  vector  of  coefficients  of  axial 
displacements. 

Second,  the  residual  is  formed: 

R  =  £(a)  +  p  (4.14) 

where  £  pertains  to  the  differential  operator  of  the  bar  equation,  that  is, 

£(u)  =  [AE(u)']' 

Third,  the  Galerkin  equations  are  formed; 

JDG(i?)d8  =  0  (4.15) 

Substituting  equation  (4.13)  and  the  bar  equation  (1.2)  into  equation  (4.14) 
yields: 

R  =  [aE(GTu)']'  +  P,(8)  (4.16) 

Substituting  equation  (4.16)  into  equation  (4.15)  3delds; 

Ji,G[AE(G’ru)']'ds  +  JdGp,(s)  ds  =  0  (4.17) 
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(4.18) 


Performing  a  single  integration  by  parts  upon  equation  (4.17)  3delds: 

AEG(GTuy  Ib  -  IdG'(  AE(GTuy)d8  +  JdGp,(s)  =  0 
Again,  removing  the  solution  vector  u  of  constant  coefficients  outside  of  the 
integral  yields: 

AEG(G’r)'u  -  JD^'[AE(GTy]ds  u  +  JdGp,(s)  =  0  (4.19) 

Recalling  that  the  axial  force  F  is  AEu',  we  define  the  boundary  term  vector 

U. 

y  =  AE  e(gTyu  Ig  (4.20A) 

the  system  stiffiiess  matrix  of  the  bar 

RA  =  JjjG'(AE(G'*'))d8  (4.20b) 

and  the  force  vector  associated  with  internal  loading  p,,  as 

r  =  lDGp,(8)  (4.20c) 

We  obtain  upon  substituting  equations  (4.20)  into  equation  (4.19), 

U-g^  +  F*s0  (4.21) 

Taking  all  the  internal  excitation  F*  terms  and  boundary  load  terms  U  to  the 
right-hand  side  of  equation  (4.21)  yields: 

=  F*  +  y  (4.22) 

Defining  the  load  vector  F^  as 

F^=  F«  +  y  (4.23) 

and  substituting  equation  (4.23)  into  equation  (4.22)  yields  the  linear  system: 

=  F^  (4.24) 

where 

5^  =  (4.25) 
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that  is,  the  system  axial  stiffiiess  matrix,  is  constructed  from  the  union  of 
the  elemental  axial  stiffiiess  matrices,  and  the  system  axial  force  vector 
is  formed  from  the  union  of  the  elemental  axial  force  vectors  f 

C.  THE  ELEMENTAL  STIFFNESS  MATRIX 

In  the  FEM  code,  the  global  (or  system)  Galerkin  FEM  equations  (4.12) 
and  (4.24)  are  actually  constructed  from  element  considerations  as  follows. 
First,  the  arch  is  divided  into  NEL  straight  beam-bar  elements  as  illustrated 
in  Figure  (2.1).  The  stress  analysis  program  contains  a  subroutine  which 
constructs  the  elemental  beam  and  bar  stiffness  matrices,  that  is,  the 
matrices  for  bending  and  k"  matrices  for  axial  stifi&iess.  The  bending  and 
axial  elemental  force  vectors,  f  and  f "  are  also  determined  within  this 
subroutine.  Figure  (4.1)  illustrates  the  degrees  of  freedom  in  which  these 
elemental  forces  (and  displacements)  act. 
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Bar  Element  Beam  Element 


Figure  4.1;  Beam  and  Bar  Element  Dcgreea  of  Freedom 

The  (4x4)  k*”  and  (2x2)  k*‘  matrices  of  the  form; 


•*^11 

«12 

^13 

^14 

A21 

k‘^ 

^22 

iC23 

K21 

k^i  = 

i-ai 

2Ci2 

^32 

ktf 

^34 

u-al 

^21 

■K22 

K:^2 

k!fi 

^44 

are  combined  to  form  a  single  (6x6)  stiffness  matrix,  This  is  accomplished 
by  redefining  the  beam  and  bar  degrees  of  freedom  in  the  following  manner: 

-  Redefine  the  bar  local  degrees  of  freedom  1'  and  2',  which  refer  to  the 
axial  displacements  at  each  end,  as  1'  and  4'  respectively. 
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-  Redefine  the  beam  local  degrees  of  freedom  1',  2',  3',  and  4',  which 
refer  to  the  lateral  displacement  and  beam  slope  at  each  end  as  2',  3', 
5',  and  6'  respectively.  The  redefined  degrees  of  freedom  are 
illustrated  by  Figure  (4.2) 

-  Place  the  respective  components  of  the  beam  matrices  and  bar 
matrices  into  the  elemental  stiffiiess  matrix  where: 


0  0  0  0 

0  k^  0  k^  ktt 

0  k^  k^  0  k^  k^ 

0  0  k^  0  0 

0  kyx  k^  0  k^  k^ 

0  k^  kS  0  k^t  , 


(4.26) 


1' 


Beam-bar  Element 


Ficnr*  4X:  Beam-Bar  Element  Degreee  oFFreedom 
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Defining  the  elemental  displacements  and  forces  in  a  similar  manner, 
the  elemental  displacement  vector  becomes: 

(5*0^  =  f (4.27) 

where  for  the  i^**  element: 

5j'  =  the  axial  displacement  at  node  1 
=  the  lateral  displacement  at  node  1 
=  the  beam  slope  at  node  1 
=  the  axial  displacement  at  node  2 
=  the  lateral  displacement  at  node  2 
^  =  the  beam  slope  at  node  2 

and  are  illustrated  in  Figure  (4.2). 

In  the  same  manner,  the  elemental  force  vector  is  redefined  as: 

(f  =  <f, fj  f,  f,  fj*',  i'  >  (4.28) 

where  for  the  i^**  element: 

fj  =  the  axial  force  at  node  1 
fg’'  =  the  lateral  force  at  node  1 
f3 ''  =  the  moment  at  node  1 
f^ ''  =  the  axial  force  at  node  2 
fg  =  the  lateral  force  at  node  2 
fg  =  the  moment  at  node  2 

also  illustrated  in  Figure  (4.2). 

With  these  developments,  the  elemental  system  of  equations  for  the 
beam-bar  element  becomes: 
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(4.29) 

or  simply 

=  (4.30) 

Prior  to  incorporation  into  the  global  matrix,  a  coordinate  transformation 
from  local  to  global  coordinates  is  imdertaken. 

D.  COORDINATE  TRANSFORMATION 

Were  all  the  elements  of  the  same  orientation  with  respect  to  one 
another  as  it  is  for  a  straight  beam,  a  global  system  of  equations  could  be 
directly  constructed.  For  the  arch,  however,  none  of  the  elements  share  the 
same  orientation.  This  necessitates  the  conversion  of  all  elemental 
displacements  and  forces  to  a  system  of  global  displacements  and  forces.  For 
a  reference  coordinate  system,  the  horizontal  and  vertical  axes  of  the  arch 
were  chosen.  (See  Figure  4.3)  Defining  the  angle  the  i^  element  makes  with 
the  horizontal  axis  as  ocj,  and  the  90°  complement  of  this  angle  as  the 
following  coordinate  relationship  between  local  and  global  "displacements" 
and  "forces"  exist. 
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^  cos(ap  +  ^  cos(p.) 
ss  cos(pi)  +  ^  COSCOj) 

% 

^'=  ^cosCaj)  +  igcos(Pj) 

^'  =  -8jco8(fj  +^co8(aj) 

S|'*  % 

fj  *'  S  fj  ‘  CO8(0Ci)  +  C08(Pj) 

fg  *'  s  fj  ‘  008(0^)  +  fj  ‘  C08{Pi) 

>'  =  £4  *  C08(aj)  +  £5  *  C080i) 

£5  =  £4  *  008(0^)  +  £5  *  C08(Pi) 


Vigan  4  J:  DupUcement  k  Fora  Coordinate  TranaformationB 


Defining  T '  as  the  transformation  matrix  for  the  i*^  element  which  is 
capable  of  performing  the  appropriate  coordinate  transformations  from  local 
elemental  coordinates  to  global  system  coordinates,  the  matrix  becomes: 


= 


cos  (ffj) 
-cos  (pj) 
0 
0 
0 
0 


COS(P_t)  0  0 

cos{a_j)  0  0 

0  10 
0  0  cos(aj) 

0  0  -cosCPj) 

0  0  0 


0 

0 

0 

cos  (Pj) 
cos  (a^) 
0 


0 

0 

0 

0 

0 

1  , 


(4.33) 


and  the  relation  between  local  and  global  "displacements"  is 


f 


< 


8i 

Ci 


’  cos  (o_^) 

cos  (pj) 

0 

0 

0 

0 

-cos  (Pjt) 

cos (a) 

0 

0 

0 

0 

82' 

0 

0 

1 

0 

0 

0 

0 

0 

0 

cos  (o^) 

cos  (p^> 

0 

•  * 

0 

0 

0 

-cos  (P^) 

cos  (ttj) 

0 

8r 

0 

0 

0 

0 

0 

1 

.  . 

(4.34) 

A  similar  relation  between  local  forces  f  and  global  forces  f  ’  exists.  The 
notation  of  equations  (4.31)  and  (4.32)  can  now  be  simplified  to 

5*'  =  r‘6‘  (4.35) 

f‘'  =  r‘f‘  (4.36) 


where: 


(5>)T  =  <5j\  6J,  5gS 

(f ’)'r  =  <fi‘,  fg’,  fj,  fg',  f6’> 


(4.37) 

(4.38) 
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E.  THE  ELEMENTAL  SYSTEM  OF  EQUATIONS 


Recall  from  equation  (4.30)  that  =  f  Substituting  equations 
(4.35)  and  (4.36)  into  equation  (4.30)  yields 

=  (4.39) 


Multiplying  both  sides  of  equation  (4.46)  by  the  inverse  of  the  transformation 
matrix,  yields 


which  simplifies  to 


(4.40) 


Since  f  Ms  an  orthogonal  matrix,  (f  *)-i  =  (f  *)T,  and  equation  (4.41)  can  be 
rewritten  as: 


(i:W'(i:*)5^  =  f‘  (4.41) 

yielding  the  elemental  system  of  equations  transformed  to  the 
horizontal/vertical  coordinate  system.  Now,  the  elemental  stiffiiess  matrices 
and  force  vectors  are  ready  for  the  construction  of  the  global  stiffness  matrix 
and  force  vector.  The  (f  ’)  term  is  the  elemental  stiffiiess  matrix  in 


terms  of  the  x-  and  y-coordinates,  and  is  denoted  as  that  is: 


=  (C')V(r’) 


(4.42) 


F.  THE  GLOBAL  SYSTEM  OF  EQUATIONS 


With  the  elemental  system  of  equations  transformed  into  the  global 
(horizontal/vertical)  coordinate  system,  the  global  system  of  equations  can  be 
formed.  The  system  stiffness  matrix  5  is  the  union  of  the  local  transformed 
stiffness  matrices  for  each  element,  thus 
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where: 


K=uk* 


(4.43) 


k’  =  (r)’^ic’'(r) 

and  the  global  (or  system)  force  vector  F  is  obtained  by  constructing  the  union 
of  the  transformed  local  force  vectors  f  *,  that  is, 

F  =  uf*  (4.44) 

Then  the  global  system  of  equations  becomes: 

55  =  F  (4.45) 

Solution  of  the  above  system  stiffness  equations  yields  the  system 
"displacements".  These  horizontal,  vertical  and  rotational  degrees  of  freedom 
"displacements"  must  be  transformed  back  to  local  axial,  lateral  and 
rotational  "displacements"  in  order  to  use  the  stress  equations  based  upon 
the  beam  and  bar  equations.  First  the  global  degrees  of  freedom  5j,  Bg,  ...  ,  5„, 

where  n=3(NEL+l),  are  related  to  the  6-  (where  i=l,2,...,6)  element 
horizontal,  vertical,  and  rotational  degrees  of  freedom  for  each  element.  The 
degree  of  freedom  for  the  i^^  element  is  given  by 

5;  =  5^  (4.46) 

where  i=l,2,...,NEL,  j=l,2,3,  and  k=3(i-l)+j.  Then  the  axial,  lateral,  and 
rotational  "displacements"  for  the  i*^  element  at  node  1  are  obtained  from 
equation  (4.31)  as 

^  coB(aj)  +  ^  cosOj) 
co8(Pj)  +  ^  cos(a,) 

and  likewise  for  the  node  2  end. 
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In  this  manner,  the  stress  at  each  end  of  each  element  can  be 
determined.  Choosing  the  greater  of  the  two  stresses  as  the  governing  stress 
of  that  element,  the  optimization  analysis  can  be  conducted  for  the  entire 
structure.  In  this  way,  the  width  dimension  b^  of  each  element  for  a 

minimized  weight  structure  is  obtained. 
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V.  PROGRAM  DESCRIPTION  AND  CAPABILITIES 


Using  the  previous  developments  of  Chapters  II,  III,  and  IV,  a 
FORTRAN  77  code  was  written  for  execution  on  the  VAX  2000  workstation. 
The  program,  named  ARCH_OPT.FOR,  was  constructed  to  be  as  fully 
interactive  with  the  user  as  possible  to  eliminate  the  need  for  editing  and 
recompiling.  The  applicability  of  the  code  is  limited  by  the  assumptions  made 
in  Chapter  H.  (i.e.,  rectangular  cross-section,  linearly  elastic  material,  etc.) 
As  illustrated  in  Figure  (5.1),  execution  of  the  program  opens  and  reads  an 
input  file,  ARCH_IN.DAT,  which  contains  information  describing  the  problem 
being  investigated.  The  x-y  coordinates  of  the  end  points  of  each  element  as 
well  as  the  element  orientation  is  determined  by  the  subroutine 
GEOMETRY.  The  subroutine  OPTIMIZATION JTOOL  contains  the 
parameter  OPTDCS,  the  optimization  decision.  With  OPTDCS=l,  DOT  is 
called  and  the  weight  optimiun  structure  is  determined  using  the  provided 
width  dimension  as  the  starting  point  of  the  opti^nization  process.  The  stress 
constraint  is  adhered  to  based  upon  the  stresses  calculated  by  the  FEM 
analysis  contained  in  the  subroutine  ARCHJSTRESS.  If  no  optimization  is 
desired,  i.e.,  OPTDCS=0,  and  the  program  computes  the  stress  distribution 
based  upon  the  input  data,  treating  the  initial  geometric  parameters  as  the 
actual  design.  With  the  data  thus  provided,  the  problem  is  solved  and  an 
output  file  named  ARCH_OUT.DAT  is  created.  The  output  file  contains  the 
problem  parameters,  the  optimized  design  variables  (width  dimensions),  and 
the  value  of  the  resulting  objective  fimction,  that  is,  the  minimum  volume. 


Fignra  B.1:  ARCH^OFT  Prognm  Stnictnre 

To  better  imderstand  the  program's  capabilities,  the  data  fields 
contained  in  ARCH_INJ)Ar  need  to  be  discussed.  The  file  is  an  unformatted 
set  of  twenty-five  numbers  separated  by  commas.  This  file  must  be  of  the 
form: 
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ANGLE,  RADIUS,  YOUNG,  YIELD,  NEL,  METHOD,  IPRINT, 
DVIBG,  DVILO,  DVIUP,  H,  CLAN,  FX,  FY,  FM,  FA, 
OPTDCS,  ITERATE,  PRCSN,  BXl,  BYl,  BMl,  BX2,  BY2,  BM2 


Table  5.1  describes  each  of  these  parameters.  For  further  clarification.  Figure 
5.2  illustrates  how  the  variables  represent  the  problem  and  the  sign 
conventions  used. 

TABLE  S.1:  ARCH_IN.DAT  FIELD  PARAMETERS 


ANGLE 

RADIUS 

YOUNG 

YIELD 

NEL 

METHOD 


IPRINT 

DVIBG 


A  real  number  from  0  to  180  representing  the  angle 
subtended  by  the  arch  (in  degrees). 

A  real  number  representing  the  length  of  the  arch. 
(Dimensions  are  arbitrary,  but  they  must  remain 
consistent  for  all  inputs!) 

A  real  number  representing  the  Y)ung's  Modulus  of  the 
arch  material. 

A  real  number  representing  the  yield  strength  of  the 
material  used.  If  a  factor  of  safety  is  desired,  it  should 
finst  be  accounted  for  and  the  resultant  design  strength 
used. 

An  integer  value  from  1  to  32  which  denotes  the  number 
of  elements  the  user  wishes  to  divide  the  arch  for  FEM 
evaluation.  The  program  is  capable  of  up  to  32  elements. 

An  integer  from  1  to  2.  This  is  a  parameter  called  by  DOT 
to  allow  the  user  to  select  which  optimization  method  is  to 
utilized. 

METHODsl  Modified  Method  of  Feasible  Directions 

METnODs2  Sequential  Linear  Programming 

NOTE:  If  the  problem  is  unconstrained,  the  BFGS 
algorithm  will  be  used  by  default  [Ref.  5,  p.  2-5] 

An  integer  from  0  to  5  used  by  DOT  to  control  the  output 
data  fi:t)m  the  DOT  optimization.  See  Appendix  C  for  the 
specific  outputs 

A  real  nxunber  which  represents  the  design  variable  1 
(width  dimension)  best  guess.  It  initializes  all  element 
width  dimension  to  the  best  guess  value.  This  establishes 
the  optimization  starting  point. 
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DVILO 

DVIUP 

H 

CLAN 


FX 

FY 

FM 

FA 

OPTIXJS 


ITERATE 


A  real  number  which  represents  desim  variable  I's  lower 
limit.  (The  lower  side  constraint  for  tne  width  dimension) 

A  real  number  which  represents  design  variable  I's  upper 
limit.  (The  upper  side  constraint  for  Sie  width  dimension) 

A  real  number  which  represents  the  constant  height 
(depth)  of  the  arch. 

An  integer  which  represents  the  number  of  the  node  at 
which  a  concentrated  load  is  to  be  applied.  This  number 
must  be  from  1  to  NEL-t-1.  If  no  concentrated  load  is 
desired,  FX,  FY,  and  FM  should  be  made  to  equal  zero. 

A  real  number  which  represent  the  magnitude  of  a 
concentrated  load  in  the  horizontal  direction.  FX  is 
applied  at  node  "CLAN”. 

A  real  number  which  represent  the  magnitude  of  a 
concentrated  load  in  the  vertical  direction.  FY  is  applied 
at  node  "CLAN". 

A  real  number  which  represent  the  magnitude  of  a 
concentrated  moment.  FM  is  applied  at  node  "CLAN*. 

A  real  number  which  represents  the  magnitude  of  a 
uniformlv  distributed  lateral  load  which  spans  the  entire 
length  of  the  arch. 

An  integer  value  which  represents  the  optimization 
"decision  such  that: 

OPTDCSsl  Optimize  the  dimension  of  the  problem 

OPTDCSs2  Do  not  optimize  the  problem.  This  choice 
will  calodate  the  stress  distribution  of  the 
arch  based  upon  the  current  problem 
dimensions,  assuming  the  width  dimension 
to  be  constant  and  equal  to  DVIBG 

An  integer  value  which  represents  the  number  of  times 
the  resisting  "optimized"  variables  are  to  be  re-entered 
into  DOT  and  the  optimization  performed  again.  Tlus 
technique  was  found  to  be  most  useful  in  refining  the 
optimized  solutions. 
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PRCSN 


B_1 


B_2 


An  integer  value  from  1  to  2.  This  parameter  allows  the 
user  to  solve  the  FEM  linear  system  of  equations  in  : 

PRCSNbI:  single  precision 

PRCSNs2:  double  precision 

An  integer  value  represent  the  boundary  condition  of  the 
arch's  first  nodal  point  such  that  if  the  value  is  0,  The  end 
is  free  to  move  in  that  degree  of  freedom.  If  the  value  is  1, 
the  end  is  fixed  in  that  degree  of  freedom. 

BXl  =  horizontal  displacement  at  point  A 
BYl  =  vertical  Replacement  at  point  A 
BMl  =  slope  of  the  beam  at  point  A 

An  integer  value  represent  the  boimdary  condition  of  the 
arch's  last  nodal  point  such  that  if  the  vRue  is  0,  The  end 
is  free  to  move  in  that  degree  of  freedom.  If  the  value  is  1, 
the  end  is  fixed  in  that  degree  of  freedom. 

BX2  =  horizontal  displacement  at  point  C 
BY2  =  vertical  displacement  at  point  C 
BM2  =  slope  of  the  beam  at  point  C 


In  summary,  for  a  given  geometry,  loading,  and  set  of  boundary 
conditions,  the  program  is  able  to  determine  the  optimum  width  dimension  of 
each  element  throughout  the  length  of  the  arch.  This  results  in  an  arch  of 
minimum  volume  (weight),  capable  of  supporting  the  given  loading.  If 
desired,  the  optimization  process  can  be  bypassed  completely.  This  results  in 
the  determination  of  the  arch  stress  distribution  based  upon  the  input 
parameters,  that  is,  the  stresses  associated  with  the  initial  design 
dimensions.  These  factors  combine  to  make  ARCHjOPT  a  useful  tool  in 
evaluating  a  variety  of  arches  in  engineering  applications. 
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Figora  &2;  ARCH_IN.DAT  Variable  Implementation 
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VI.  PROGRAM  VERIFICAriON 


To  verify  the  finite  element  code  used  in  the  optimization  investigation, 
several  non-optimum  problems  of  beam  and  arch  structures  with  known 
analytical  solutions  were  solved  iising  the  code.  The  code  solution  of  these 
problems  also  served  the  purpose  of  establishing  a  relation  between  the 
number  of  elements  used  in  the  model  and  the  accuracy  of  the  method.  With 
a  "yardstick"  thus  provided,  "ARCH_OPT.FOR"'s  capabilities  for  accurate 
modeling  of  beams  and  arches  was  assumed. 

The  first  verification  problem  was  a  cantilever  beam  subjected  to  a 
concentrated  end  load,  illustrated  by  Figure  (6.1).  Gere  and  Timoshenko 
[Ref.  6,  p.  737]  give  general  formulas  for  the  lateral  displacements  and  beam 
slopes  for  this  case,  as: 

V  =  Px2(3L-x)/6EI  (6.1a) 

v' =  Px(21-xy2EI  (6.1b) 

Using  the  parameters: 

P  =  1000  Ibf  h  =  3.0  inches 

L  =  45  inches  £  =  30  x  10®  psi 

b  =  1.5  inches  I  =  bh®/12 

the  lateral  displacements  and  beam  slopes  at  the  midpoint  and  free  end  were 
calculated.  ARCHjOPT  was  nm  for  this  beam  structure  using  an  angle  of 
45.0  X  10'®  radians  and  a  radius  of  10®  inches  to  approximate  the  straight  45 
inch  beam  length.  For  a  four  element  approximation.  Table  (6.1)  compares 
the  FEM  solution  to  the  analytic  solution. 
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Node 

Fixed  End 

Mid  Point 

Free  End  1 

Analytic  5 
FEM  5 

O.OOOE+00 

O.OOOE+00 

9.375E-02 

9.375E-02 

3.000E-01 

3.000E-01 

%  Error 

fixed 

O.OOE+00 

O.OOE+00 

Analytic  5' 
FEM  5’ 

O.OOOE+00 

O.OOOE+00 

7.500E-03 

7.499E-03 

l.OOOE-02 

9.999E-03 

%  Error 

fixed 

1.33E-02 

l.OOE-02 

Max  %  Error  1+33E-02:S 

where  percent  error  is  defined  as, 

%  Error  *  100  x  (5.^^  -  5pE^)/  (6.2) 


The  stress  corresponding  to  the  same  points  of  interest  was  calculated 
using  equation  (3.1).  Recall 


=  Mc/I 
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which  in  terms  of  the  beam  equation  at  the  i^**  node  becomes 

(a  ).  ^  (EIv"h/2)/(b;h3/12) 


or  simply 


(a„).  =  6EIv"/bih2 


(6.3) 


From  equation  (4.1)  of  the  FEM  development,  v  =  Q’^y,  hence 

v"  =  (QTy)" 

V"(x^)  =  (Q'^rv  =  (-6/l2)vj  +  (4/1)y/  (6.4) 

Substituting  equation  (6.4)  into  equation  (6.3)  yields 

(o„)i  =  (6EL^jh2)[-6vj/l2  +  4v//l]  (6.5) 

where  Vj  and  Vj'  are  the  lateral  displacement  and  slope  at  and  are  obtained 
from  5^',  8^',  8^'  and  8^'  of  equation  (4.31).  The  stresses  of  equation  (3.1)  were 
then  compared  to  those  calculated  by  the  ARCHjOPT  using  equation  (6.5). 
The  results,  summarized  in  Table  (6.2),  show  a  maximum  difference  of 
5.00x10  '*%  between  the  exact  and  FEM  solution. 


TABLE  6JJ;  VERinCATION  PROBLEM  #1  SUMMARY  OF  STRESSES 

Verification  Problem  #1 
(Stresses) 


Node 

Fixed  End 

Mid  Point 

Free  End  1 

Analytic  <t 
FEM  or 

2.000E+04 

1.999E+04 

l.OOOE+04 

9.999E+03 

O.OOOE+00  j 
9.712E-06 

%  Error 

5.00E-04 

5.00E-04 

N/A  \ 

Max  %  Error 


The  second  verification  problem  considered  was  that  of  a  prismatic  bar, 
fixed  at  one  end  and  subjected  to  an  axial  concentrated  load  at  the  free  end  as 
shown  in  Figure  (6.2). 
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Fiffoiv  fti2:  Verification  Problem  #2 

The  parameters  used  were: 

P  =  1000  Ibf  h  =  3.0  inches 

L  =  45  inches  E  =  30  x  10®  ps> 

b  =  1.5  inches  I  =  bh®/12 

From  the  Bar  Equation  (1.2),  we  hs  ve: 

(AEu0'  =  p,(x) 

AEu'  =  F(x)  (6.6) 

and  the  normal  stress  developed  in  a  bar  due  to  axial  loading  is: 

=  F(x)/A  (6.7) 

Substituting  equation  (6.6)  into  equation  (6.7)  yields 

o„  =  AEu7A 

or  simply 

On  =  Eu'  (6.8) 

Substituting  equation  (4.18)  into  equation  (6.8)  yields  the  FEM  solution, 

o„  =  E(GTu)'  (6.9) 
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From  the  FEM  development, 


(GV  =  (Ui/l) 

which  when  substituted  into  equation  (6.9)  yields  the  stress  at  the  i‘^  node, 

(o„)i  =  Eui/l  (6.10) 

where  Uj,  the  axial  displacement  at  Xj,  is  obtained  from  6['  and  5^  of  equation 

(4.31).  Again  using  the  midpoint  and  free  end  as  the  locations  for  conducting 
the  displacement  and  stress  analysis,  the  FEM  code  results  are  compared  to 
those  generated  by  the  analytic  solutions.  These  results  are  listed  in  Table 
(6.3).  They  show  no  difference  between  the  FEM  approximation  and  the 
analytic  solution. 

TABLE  « VERIFICATION  PROBLEM  #2  SUMMARY  OF  RESULTS 


Verification  Problem  #2 
(Displacements) 


1 _ Node _ 1 

Fixed  End 

Mid  Point 

Free  End  I 

Analytic  5 

O.OOOE+00 

1.667E-04 

3.333E-04 

FEM  5 

O.OOOE+00 

1.667E-04 

3.333E-04 

%  Error 

fixed 

O.OOE+00 

O.OOE+00 

Analytic  5' 

O.OOOE+00 

2.222E+02 

2.222E+02 

FEM  5' 

O.OOOE+00 

2.222E+02 

2.222E+02 

%  Error 

fixed 

O.OOE+00 

O.OOE+00  1 

Max  %  Error  0,015^1^00:1; 


The  third  verification  problem  chosen  was  a  cantilever  beam  with  a 
concentiated  moment  at  the  free  end,  illustrated  by  Figure  (6.3).  Again,  from 
Gere  and  Timoshenko  [Ref.  6,  p.737],  the  analytic  solution  for  the 
displacements  and  slopes  of  this  particular  problem  are  given  by  the 
equations: 


V  =  M„x2/2EI 


(6.11a) 


Figure  ftJ:  Verification  Problem  #3 


U.^ing  the  parameters: 

=  10,000  Ibf 
L  s  45  inches 
b  =  1.5  inches 


h  =  3.0  inches 

E  =  30  X  10®  psi 
I  =  bh3/12 


the  displacements  and  slopes  were  determined  for  the  points  of  interest. 
These  results  are  compared  to  the  4-element  FEM  approximations  in  Table 
(6.4). 


TABLE  6.4:  VEROTCATION  PROBLEM  #3  SUMMARY  OF  DISPLACEMENTS 


Verification  Problem  j({3 


Node 

Fixed  End 

Mid  Point 

Free  End 

Analytic  5 

O.OOOE+00 

2.500E-02 

l.OOOOE-01 

FEM  5 

O.OOOE+00 

2.500E-02 

9.9999E-02 

%  Error 

fixed 

O.OOE+00 

l.OOE-03 

Analytic  5' 

O.OOOE+00 

2.222E-03 

4.444E-03 

FEM  5’ 

O.OOOE+00 

2.222E-03 

4.444E-03 

%  Error 

fixed 

O.OOE+00 

O.OOE+00 

Max  %  Error 


l.OOE-03 


A  comparison  of  the  FEM  approximations  and  analytical  solutions 
again  show  virtually  no  disparity.  These  results  are  presented  in  Table  (6.5). 
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TABLE  «JS:  VERIFICATION  PROBLEM  #3  SUMMARY  OF  STRESSES 
Verification  Problem 


(Stresses) 

Node 

Fixed  End 

Mid  Point 

Free  End 

Analytic  <r 
FEM  <T 

4.444E+03 

4.444E+03 

4.444E+03 

4.444E+03 

4.444E+03 

4.444E+03 

%  Error 

O.OOE+00 

O.OOE+00 

O.OOE+00 

Max  %  Error 


The  final  verification  problem  is  based  upon  an  example  arch  problem 
presented  by  Gere  and  Timoshenko  [Ref.  6,  p.616].  They  demonstrate  how 
the  umt-load  method  can  be  used  to  calculate  the  horizontal  displacements  of 
the  problem  illustrated  in  Figure  (6.4).  The  following  formula  was  obtained; 

5h  =  PR3/2EI  (6.12) 


Figure  6.4:  Arch  Verification  Problem  *4 


For  the  parameters: 


P  =  1000  Ibf  e  =  90“ 

R  =  45  inches  E  s  30  x  10^  psi 

b  =  1.5  inches  I  =  bh^/12 

h  =  3.0  inches 

the  horizontal  deflection  is  found  to  be  0.4500  inches.  The  horizontal 
deflection  from  the  4-element  FEM  approximation  is  0.4470  inches,  an  error 
of  0.66%. 

In  all  of  the  verification  problems,  the  percent  differences  between  the 
values  obtained  by  the  approximate  FEM  method  and  the  exact  solutions 
were  in  all  cases  less  than  0.66%.  Satisfied  that  the  program  was  producing 
very  good  data,  the  investigation  to  obtain  optimum  structures  was  pursued. 
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VII  CASE  STUDIES 


Recognizing  an  arch  to  be  but  another  structural  means  of  transferring 

a  load  from  one  point  to  another,  the  desire  to  compare  the  efficiency  of  the 

optimized  arch  to  that  of  a  traditional  structure  drove  the  first  two  cases 

studied.  Given  the  problem  of  transferring  the  load  at  B  to  A,  as  illustrated 

in  Figure  (7.1a),  niimerous  structures  could  be  used.  For  brevity,  only  the 

frame  (7.1b)  and  arch  (7.1c)  will  be  studied.  Given  the  parameters: 

E  =  30x10^  psi  h  =  2  inches 

Sy  =  52,000  psi  a  =  32  inches 

I  =  bh^/12  b  =  32  inches 

only  the  width  dimension  for  each  case  will  be  allowed  to  vary.  In  this  way,  a 
volume  comparison  of  each  structure,  hence  a  measure  of  the  relative 
efficiency  of  the  structure,  may  be  made. 

For  the  non-optimized  frame  shown  in  Figure  (7.1b),  the  base  (or  width) 
dimension  is  considered  constant  throughout.  In  order  to  keep  the  maximum 
stress  below  the  yield  strength  of  the  arch  material,  we  have 

or 

M„„c/I  ^  52,000  psi  (7.1) 

The  maximum  moment  occurs  uniformly  along  the  vertical  member  of  the 
frame,  hence 

[(2000  lbf)(32  inX2  in/2)Mb(2  in)3/12]  ^  52,000  psi 
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or  simply 


(7.2) 


1.846  in.  sb 

The  total  volume  of  the  frame  is 

Volume  s  (32  inX2  inXl.846  in)  -f  (32  inX2  inXl.846  in) 
or 

Volume  *236.3  in®  (7.3) 

This  volume  will  be  the  basis  upon  which  the  following  optimized 
volumes/weights  hence  efficiencies  will  be  based. 


(b)  (c) 

PiforB  74:  Methods  of  l^mnsforring  a  Load 

A.  CASE  1:  THE  OPTIMIZED  FRAME 

Given  the  problem  presented  in  Figure  (7.1),  an  optimized  frame  of 
equal  load  bearing  capabilities  was  sotight.  First,  the  frame  was  divided  into 
its  vertical  and  horizontal  members.  The  vertical  member  is  subjected  to  the 
concentrated  moment,  aP,  at  C,  resulting  in  a  uniform  moment  and  hence 
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uniform  stress  along  the  member.  Consequently  the  vertical  member  has  a 
uniform  width  dimension  of  1.846  inches  as  previously  determined. 

The  horizontal  member  is  a  cantilever  beam  subjected  to  a 
concentrated  end  load.  Since  the  moment  along  BC  varies  linearly  from  0  at 
B  to  aP  at  C,  for  a  constant  g^,,=Sy,  the  width  must  also  vary  linearly.  Thus 

the  width  varies  linearly  fmm  0  inches  at  the  right  end  to  1.846  inches  at  C. 
The  volume  for  this  frame  is: 

Volume  =  (32  inX2  inXl.848  in)  +  (.5X32  inX2  inXl.848  in) 

which  is; 

Volume  =  177.4  in®  (7.4) 

The  volume  of  the  horizontal  member,  BC,  was  then  optimized  using  a 
4-element,  8-element,  and  12-element  discretization.  Table  (7.1)  illustrates 
the  optimized  volumes  of  each  of  these  solutions.  The  percent  difference 
between  the  4-element  and  8-element  solution  was  found  to  be  less  than 
11.5%  and  that  for  the  8-element  and  12-element  solution  to  be  less  than  3%. 
In  the  interest  of  solving  many  cases,  it  was  decided  to  solve  all  future  case 
studies  with  a  12-element  discretization,  treating  the  12-element  model  as 
producing  grid  independent  resiilts. 

The  optimized  results  for  the  horizontal  member,  given  in  Table  (7.1) 
show  a  total  member  volume  of  64.63  in®.  The  volume  of  the  vertical  member 
remains  the  same  as  for  the  non-optimized  structure,  hence  the  total  volume 
of  the  optimized  frame  is 

Volume  =  (32  inX2  inX  1.848  in)  +  64.63  in® 

=  182.9  in®  (7.5) 
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This  represents  22.6%  less  volume,  hence  22.6%  less  weight  than  the  non- 
optimized  frame. 


TABLE  7.1:  CASE  1  SUMMARY  OF  RESULTS 


Element 

Height 

hnchesi 

Length 

hnchesi 

Base 

hnchesi 

Volume 
{cubic  in.i 

Node 

Stress 

{nsii 

2.000E+00 

2.687E+00 

l.B4BE-t-00 

9.B56E4O0 

5.192-1-04 

2000E44)0 

2.667E+00 

1.6g6E400 

B.040E-H)0 

5.192-104 

2.000E+00 

2.667E+00 

1.543E+00 

8.22gE-t-00 

5.182-1-04 

4 

2  000E-h00 

2.667E+00 

1  38gE400 

7  40BE+00 

5.182-104 

5 

2.000E+00 

2.667E44}0 

1.23eE+00 

6.5622-1^0 

5.182404 

6 

2000E+00 

2.667E+00 

i.0B2E-«-00 

5.7712-1-00 

6 

5.182+04 

2.000E+00 

2.e87E+00 

e.313E-01 

4.967E-K)0 

n 

5.172+04 

2.000E+00 

2.667E+00 

7.7B1E-01 

4.1342-1-00 

mm 

6.152+04 

2.000E+00 

2.667E+00 

6.318E-01 

3.3B6E-l4}0 

5.162+04 

2  000E+00 

2.667E+00 

4.676E-01 

2  404E-1-00 

5.072+04 

2  000E-H)0 

2.667E+00 

3.196E-01 

1.705E-1-00 

5.132404 

2  000E+00 

2  667El4)0 

2.000E-01 

1.067E-1-00 

6.072+04 

12 -element 

£  Volume: 

4.002+04 

B-eleraenl  E  Volume:  B.BABE-l-Ol 
4-element  £  Volume:  7.3B5E+01 


B.  CASE  2:  THE  OPTIMIZED  CANTILEVER  ARCH 

Given  the  circumstances  and  parameters  of  Figure  (7.1),  a  cantilever 
circular  arch  (Fig.  7.1c)  was  employed  to  perform  the  same  function, 
transferring  the  given  load  at  point  B  to  point  A.  The  resulting  dimensions 
and  stresses  of  the  optimization  are  presented  in  Table  (7.2).  These  results 
illustrate  what  one  would  have  expected,  the  width  dimension  is  reduced 
until  the  local  stress  approaches  the  yield  strength  of  the  material.  The  total 
volume  of  the  arch,  128.3  in®  is  46.3%  less  than  the  non-optimized  frame  and 
29.9%  less  than  the  optimized  frame.  In  moving  structures  where  higher 
weights  mean  higher  operating  costs,  savings  such  as  these  can  become 
significant. 
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TABLE  7Je:  CASE  2  SUMMARY  OF  RESULTS 


Clement 

Base 

hnchesl 

Volume 
{cubic  in.{ 

Node 

Stress 

h)8i{ 

1 

2.000E+00 

4.1B6E-tOO 

1.B64C+00 

1.561E401 

5.20E4O4 

2 

2.000E400 

4.1B6E+00 

1.B42E400 

1.B42E-I01 

5.20E+04 

3 

2.000E400 

4.1B6E400 

1.705E400 

1.503E401 

5.20E404 

4 

2.000E+00 

4.1B6E400 

1.71BE-«O0 

1.43BE401 

B.20E404 

5 

2.000E-K)0 

4.186E400 

1.612E400 

1.350E401 

5 

5.10C404 

6 

2.000E400 

4.1BBE400 

1.477E4O0 

1.237E401 

6 

5.1BE404 

7 

2.000E400 

4.1B6E4O0 

1.325E4O0 

l.lOOE+01 

■■ 

5.16E+04 

B 

2.000E-fO0 

4.1B6E400 

1.166E400 

9.7B2E4O0 

5.16E-K)4 

0 

2.000E-K)0 

4.1B6E400 

8.316E-01 

7.7B9E4O0 

5.05E-fO4 

10 

2.000E-K)0 

4.1B6E400 

7.126E-01 

5.B6BE400 

5.19E404 

11 

2.000E400 

4.1B6E400 

5.16BE-01 

4.327E-H}0 

4.B5E404 

12 

2000E400 

4.1B6E400 

a656E-01 

3.061E+00 

12 

3.47E404 

X  Volume: 

13 

1.79E+02 

C.  CASES:  THE  OPT] 
END  LOAD 


SD  CANTILEVER  ARCH  WITH  AXIAL 


To  appreciate  how  the  stress  distribution  follows  the  moment 
distribution  of  the  member,  the  problem  of  Case  2  was  modified  keeping  the 
same  parameters  as  outlined  previously,  by  changing  the  direction  of  the  load 
as  shown  in  Figure  (7.2). 
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For  this  case,  the  moment  at  any  point  is  given  by  the  expression 

M  =  PR(l-sin(e))  (7.6) 

as  opposed  to  the  moment  distribution  of  Case  2  where  the  moment  along  the 
length  of  the  arch  is 

M  =  PRcos(e)  (7.7) 

Assuming  a  unit  moment,  that  is  PR=1,  then  from  Figure  (7.3),  one  may  see 
that  the  area  under  Case  3's  moment  curve  is  significantly  less  than  the  area 
under  Case  2's  moment  curve.  One  would  expect  this  to  correspond  with  the 
need  for  less  material  due  to  less  applied  force  (in  this  case,  moment). 
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Moment  Diagram 


The  resulting  optimized  arch  of  this  case  has  the  dimensions  outlined 
in  Table  (7.3)  and  a  noticeably  smaller  volume  of  81.39  in®  than  the  128.3  in® 
of  Case  2.  In  fact,  the  reduction  in  volume  is  the  same  as  the  reduction  in 
areas  of  the  moment  diagrams.  This  result  collaborates  with  the  previous 
observation.  Hence  one  may  see  that  the  normal  stress  is  predominantly  due 
to  bending  and  essentially  follows  the  moment  distribution  of  the  structure  in 
question. 
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TABLE  7A:  CASE  3  SUMMARY  OF  RESULTS 


EUement 

■Una  SB 

Node 

Stress 

h>8i{ 

1 

2.000E400 

4.1B6E-fOO 

1.846E4O0 

1.545E401 

1 

5.20E-«O4 

2 

B.OOOEi-OO 

4.1B6E+00 

1.6B7E400 

1.412E401 

2 

5.19E404 

3 

2.000E400 

4.1BBE400 

1.435E4O0 

1.201E401 

3 

5.20E4O4 

4 

2.000E400 

4.1BBE400 

1.196E400 

l.OOlE+01 

4 

6.1BE404 

5 

2.000E-^0 

4.1BBE400 

0.703E-01 

B.123E400 

5 

5.1BE404 

6 

2.000E400 

4.1B6E-tOO 

7.674E-01 

B.341E400 

6 

6.1BE404 

7 

2.000E400 

4.1BSE400 

5.676E-01 

4.752E4O0 

7 

5.19E404 

B 

2.000E400 

4.1B6E400 

4.00gE-01 

3.3BBE+00 

B 

5.1BE404 

0 

2.000E-fO0 

4.1B6E400 

2.611E-01 

2.1B6E4O0 

B 

5.16E404 

10 

2.000E-^0 

4.1B6E400 

2.000E-01 

1.B74E4O0 

10 

4.7SE404 

11 

2.000E400 

4.1BbE400 

2.000E-01 

1.674E4O0 

11 

2.13£-(04 

12 

2000E-M)0 

4.1B6E400 

2.000E-01 

1.674E+00 

12 

9.10E4C3 

£  Volume; 

13 

5.00E+03 

D.  CASE  4:  THE  CANTILEVER  ARCH  UNDER  A  DISTRIBUTED 
LOAD 


This  case  involves  applying  a  uniformly  distributed  load  acting  radially 
inward  on  a  cantilever  arched  segment  as  pictured  in  Figure  C7.4) 


Picnra  tA:  C*m  4  Problvm  Gaometry 
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where: 

E  =  30x10®  psi  h  =  2  inches 

Sy  =  52,000  psi  R  =  32  inches 

I  =  bh®/12  0  =  90  degrees 

The  results  of  the  optimization  contained  in  Table  (7.4)  illustrate  how 
the  critical  constraints  shift  from  the  maximum  allowable  stress  to  the 
minimum  allowable  width  dimension.  Even  though  the  moment  in  the 
structure  has  diminished  and  the  stress  no  longer  approaches  the  3deld 
strength  of  the  material,  the  geometric  constraint  prevents  the  cross-section 
from  becoming  so  thin  that  the  bar/beam  assumption  is  no  longer  valid.  The 
variation  in  the  width  dimensions  appears  to  be  almost  logarithmic  alluding 
to  the  complexities  involved  with  arched  segments  under  uniformly 
distributed  loads.  Had  the  arch  not  been  optimized  and  assuming  the 
maximum  width  dimension  of  Table  (7.4)  to  represent  the  uniform  width 
dimension  of  the  non-optimum  arch,  such  that 

Volume  =  [(R^)2  -  (Rj)2](e  radiansXb^„)/2  (7.8) 

then  the  total  volume  of  the  non-optimum  arch  would  have  been 

Volume  =  209.1  in.® 

Hence,  the  optimized  volume  of  101.6  in®  is  51%  smaller  than  that  of  the 
non-optimized  arch. 
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TABLE  7.4;  CASE  4  SUMMARY  OF  RESULTS 


Element 

Height 

hnchesl 

Node 

Stress 

h)8i{ 

1 

2.000E+00 

4.1B6E+00 

2.0BOE400 

1.741E+01 

5.21E-K04 

2 

2.000E-K)0 

4.1B6E-f00 

1.930E400 

1.818E401 

5.20E-f04 

3 

2.000E400 

4.1B6E4O0 

1.782E+O0 

1.475E401 

5.20E-fO4 

4 

2.000E400 

4.1B6E400 

1.552E4O0 

1.2BQE+01 

8.20E4O4 

5 

2.000E400 

4.1B6E-fOO 

1.315E4O0 

l.lOlE+01 

5 

5.20E-K)4 

8 

2.000E400 

4.1B6E-H}0 

1.087E4O0 

B.933E+O0 

5.20E4O4 

E.OOOE+00 

4.1B6E+00 

6E13E-01 

B.B78E+00 

5.20E-f04 

2.000E-tO0 

4.1B6E400 

&.912E-01 

980E400 

5.19E404 

2.000E400 

4.1B6E400 

3.001E-01 

3.266E-K)0 

9 

5.18E404 

2.000E400 

4.1B6E+O0 

2.21BE-01 

1.B87E+O0 

10 

8.1BE-K)4 

2.000E400 

4.1B6E400 

2.000E-01 

1.674E4O0 

11 

2.52E404 

2.000E-fO0 

4.1B8E4O0 

2.000E-01 

1 674E+00 

12 

B.82E-K)3 

£  Volume: 

13 

l.lOE+03 

E.  CASE  6:  THE  SIMPIY  SUPPORTED  ARCH 


This  case  involved  the  optimization  of  a  simply  supported  arch 
subjected  to  a  lateral  load  at  the  midpoint,  illustrated  in  Figure  (7.5).  The 
arch  is  a  first  order  statically  indeterminate  structure  subject  to  the  following 
parameters: 

E  s  30x10^  psi  hs  2  inches 

S  =  52,000  psi  R  s  32  inches 

I  s  bh^/12  6  s  180  degrees 
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F>  =  8000  lt>f^ 


In  order  to  take  advantage  of  symmetry,  the  problem  was  divided  along  the 
axis  of  symmetry  with  the  following  boundary  conditions  imposed  on  the 
symmetry  end, 

(EIv")'  =  P/2 
EIv'  =  0 
u  =  0 

The  results  of  the  optimization  are  summarized  in  Table  (7.5).  Here 
again,  the  weight  savings  of  the  optimized  arch  of  263.2  in^  over  that  of  the 
non-optimized  arch,  as  defined  by  equation  (7.8),  of  628.1  in^  is  approximately 
58%. 


-57- 


TABLE  7^:  CASE  6  SUMMARY  OF  RESULTS 


Element 

HS5<I  SH 

■InSiB 

Node 

Stress 

hail 

1 

2.000E400 

4.1B6E-K)0 

S.130E-01 

4.2B5E4O0 

B.0BE-tO3 

2 

2.000E400 

4.186E4O0 

B.060E-01 

7.5g2E400 

6B0E-K)4 

3 

2.000E400 

4.1B6E-K)0 

i.ll7E-fO0 

B.351E-K)0 

5.00E+04 

4 

2.000E4U0 

4.1BBE400 

1.47aE400 

1.23BE-f01 

5.20E4O4 

5 

2.000E400 

4.1BBE400 

l.lOBE-fOO 

1.001E-K)1 

5 

5.20E404 

6 

2.000E400 

4.1B6E^  30 

1.127E4O0 

0.435E400 

6 

5.20E404 

2.000E400 

4.1B6E+00 

1.531E4O0 

1.2B2E+01 

4.22E404 

2.000E-M)0 

4.1B6E400 

6.785E-01 

4.B26E4O0 

4.01E404 

2.000E-KO0 

4.1B6E-K)0 

5.03SE-01 

4.069E400 

4.93E403 

2.000E400 

4.1BBE400 

1.346E4O0 

1.127E-f01 

10 

5.20E4O4 

2.O00E-fO0 

4.1B6E4O0 

2.206E400 

1.B47E401 

11 

5.20E4O4 

12 

2000E400 

4.1B6E400 

3.124E400 

2.615E-K)1 

12 

B.20E4U4 

X  Volume: 

13 

S.20E+O4 

F.  CASE  6:  THE  FIXED-FIXED  ARCH 

To  determine  the  effect  of  additional  redundancy,  the  next  case 
involved  the  optimization  of  a  fixed-fixed  arch  subjected  to  a  lateral  load  at 
the  midpoint,  illustrated  in  Figure  (7.6)  and  subject  to  the  following 
parameters: 

E  =  30x10^  psi  h  s  2  inches 

Sy  =  52,000  psi  R  =  32  inches 

I  =  bh^/12  0  =  180  degrees 
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The  fixed  supports  add  two  additional  redundant  moments  at  the  supports 
when  compared  to  the  previous  simply-supported,  simply-supported  arch  of 
case  5.  Again  taking  advantage  of  symmetry,  the  problem  was  divided  along 
the  axis  of  symmetry  with  the  following  boundary  conditions  imposed  on  the 
symmetry  end, 

(EIv")'  =  P/2 
EIv'  =  0 
u  =  0 

The  resulting  optimized  arch  has  a  total  volume  of  210.3  in^  as 
summarized  in  Table  (7.6).  The  non-optimum  arch  has  a  volume  of  533.0  in^ 
assuming  a  constant  arch  width  corresponding  to  the  optimized  arch's 
maximum  width  of  2.651  inches.  Here,  the  weight  savings  of  the  optimized 
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arch  over  that  of  the  non-optimized  arch  is  approximately  60%.  It  is  also 
worth  particular  note  that  though  the  loading  and  geometry  of  Cases  6  and  6 
are  the  same,  differing  only  with  regards  to  the  boimdaiy  conditions  imposed, 
Case  6  is  20%  lighter.  This  is  because  a  fixed-fixed  structure  is  more 
statically  indeterminate  than  the  simply  supported  structure.  Hence  we  see 
here  that  the  more  redundant  structure  is  the  more  efficient  member.  This 
illustrates  one  of  the  reasons  why  fixed-fixed  structures  are  preferred  in 
construction. 

TABLE  7.6;  CASE  6  SUMMARY  OF  RESULTS 


Element 

Height 

Unchesi 

Length 

Uncheal 

Base 

ilncheal 

Volume 
{cubic  ln.( 

Node 

Streas 

_ i£Sil _ 

1 

2000E+00 

4.186E+00 

1.727E+00 

1.44BE-i-01 

5.20E+04 

2 

2.000E+00 

4.186E+00 

g.458E-01 

7.91BE+00 

5.00E+04 

3 

2  000E+00 

4.186E-l'00 

2  710E-O1 

2.276E+00 

4.62E+04 

4 

2.000E4-00 

4.186E400 

eS04E-01 

5.445E+00 

B.20E+04 

5 

2.000E-»-00 

4.186E+00 

1.001E400 

9.133E+00 

5 

4.94E+04 

6 

2000E+00 

4.186E+O0 

7.899E-01 

B.613E4P0 

6 

5.20E+04 

2.000E+00 

4.186E4-00 

1.062B+00 

B.891B4O0 

mm 

5.1BE+04 

2  000E+00 

4.1BeE-f00 

7.539E-01 

6.311E+00 

4.21E+04 

2.000E-K)0 

4.1B6E+00 

2.585E-01 

2.164E-H)0 

5.20E+04 

2,000E+00 

4.1B6E-t^0 

9  328E-01 

7.B09E+00 

5.13E+04 

2.000E+00 

4.1B6E+00 

I.753E4O0 

1.46BE+01 

5.24E+04 

2.000E+00 

4.186E+00 

2.e51E4O0 

2.219E4ni 

5 19E+04 

£  Volume; 

5.I9E+04 

VIII  CONCLUSIONS 


The  conclusions  of  this  study  are: 

-  The  stress  analysis  based  upon  the  bar/beam  model  yielded  good 
results  with  percent  deviations  from  known  anal3rtic  solutions 
ranging  between  0.1  and  1.5%.  Hence,  the  bar/beam  element  model 
is  a  viable  technique  in  the  approximation  of  arch  structures. 

-  The  DOT  optimization  software  was  able  to  utilize  the  bar/beam 
modeled  stress  analysis  to  efficiently  determine  weight  optimum 
arch  structures. 

-  The  optimization  demonstrates  how  structures  which  are  more 
statically  indeterminate  (redundant)  are  likewise  more  efficient  than 
identical  structures  under  identical  loading. 

-  The  weight  optimization  of  a  structure  is  available  and  effective  for 
all  types  of  problem  boundary  conditions. 

It  should  be  noted  once  again  that  this  is  an  initial  investigation  into 
the  weight  optimization  of  arches.  Numerous  opport\mities  exist  for  the 
expansion  of  the  basic  assumptions  made  in  this  study.  This  investigation 
only  considered  the  optimization  of  structures  of  the  form: 

I  =  kA 
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The  general  form  of  this  type  of  optimization  is  such  that 

I  =  kA« 

where  n  =  1,  2,  or  3  depending  upon  which  cross-sectional  dimension(s)  is 
defined  as  the  design  variable.  Some  of  the  possibilities  for  future  research 
include: 


-  Allowing  only  the  cross-sectional  height  dimension,  h,  vary  while 
holding  all  other  parameters  constant.  (n=3) 

-  Allowing  both  the  height  and  width  dimensions  to  vary 
proportionally,  while  holding  all  other  parameters  constant.  (n=2) 

-  Allowing  the  radius  of  curvature  of  the  arch  (its  center-line  shape)  to 
vary. 

-  Optimizing  the  arch  using  engineering  cross-sections  such  as  box 
beams,  I-beams,  circular  cross-section,  etc. 

-  Incorporating  additional  constraints  (such  as  arch  maximum  height 
limitations,  buckling  constraints,  crippling  constraints,  etc,)  in  order 
to  expand  the  model;  thereby  enabling  the  model  to  solve  a  greater 
variety  of  problems. 
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APPENDIX  A 


JUSTIFICATION  FOR  OMITTING  SHEAR  STRESSES 

The  shear  stress  distribution  through  a  beam  of  rectangular  cross- 
section  has  a  parabolic  distribution  along  the  height  of  the  member.  The 
maximum  shear  stress,  located  at  the  neutral  axis  of  the  beam,  is 

W  =  1-5V/A  (A.1) 

where  is  the  maximum  shear  stress,  V  is  the  shear  force,  and  A  is  the 
cross-sectional  area  of  the  beam.  [Ref.  6,  p.  229] 

The  normal  stress  due  to  bending  is  given  by  the  equation 

<T^  =  Mc/I  (A.2) 

where  is  the  maximum  normal  stress,  M  is  the  bending  moment,  and  I  is 
the  cross-sectional  moment  of  inertia  which  for  this  case  is  bh^/12  where  b 
and  h  are  the  width  and  height  respectively  of  the  cross-section. 

Redefining  the  normal  stress  in  terms  of  the  cross-sectional  dimensions 

yields 

=  M(h/2)/(bh3/12) 
or 

=  6M/hA  (A.3) 

The  ratio  of  the  maximum  shear  stress  to  the  normal  stress  due  to 
bending,  is  denoted  by  r  and  given  by  the  expression: 

r  =  (A.4) 
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Substituting  equations  (A.1)  and  (A.3)  into  equation  (A.4)  yields 

r  =  (1.5V/Ay(6M/hA) 


or 


r  =  Vh/4M  (A.5) 

For  the  cases  investigated  in  this  study,  the  maximum  value  r  can  attain  is 
when  the  loading  is  that  of  a  uniformly  distributed  load,  p^.  Then,  where: 

V  =  PyL  (A.6) 

M  =  PyL2/2  (A.7) 

which  upon  substitution  into  equation  (A.8)  yields 

r  =  (PyI.)h/4(pyL2/2) 

which  simplifies  to 

r  =  h/2L  (A.8) 

The  use  of  the  beam  equation  requires  the  length  of  the  beam  to  be  at  a 
minimum  ten  times  the  height,  that  is: 


L^lOh 


(A.9) 


To  maximize  the  value  of  r,  let  L  equal  lOh,  the  minimum  allowable  length. 
Cabstituting  this  value  of  L  into  equation  (A.8)  yields 

r^h/2(10h) 


or  simply 


r^l/20 


(A.10) 


Hence,  the  muTimuTn  shear  stress  accounts  for  less  than  five  percent  of  the 
bending  stress  developed  in  the  structure.  Five  percent  is  high  considering 
this  analysis  over-assumed  the  value  of  the  shear  stress  by  assigning  the 
maximum  shear  stress  to  the  entire  cross-section  of  the  beam.  Moreover,  at 
the  outermost  fibers  where  is  a  maximum,  the  shear  stress  is  zero. 
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Therefore,  under  the  circumstances  of  this  study,  the  addition  of  shear 
stresses  was  deemed  to  be  unwarranted. 
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APPENDIX  B 


DOT  USERS  MANUEL,  SECTION  2.1 

DOT  WITH  APPLICATION  PROGRAMS _  _ 


2.1  CALLING  STATEMENT _ 

DOT  is  invoked  by  the  following  FORTRAN  calling  statement 
in  the  user’s  program; 

CALL  DOT  ( INFO,  METHOD,  IPRINT,  NDV,  NCON,  X, 

•  XL,  XU.  OBJ,  MINMAX,  G,  RPRM,  IPRM,  WK,  NRWK, 
•IWK,NRIWK) 
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APPENDIX  C 


DOT  USERS  MANUEL,  SECTION  2^ 


2.2  PARAMETERS  IN  THE  CALLING 
STATEMENT 


T&ble  2-1  lists  the  parameters  in  the  calling  statement  to  DOT. 
Where  arrays  are  defined,  the  required  dimension  size  is  given 
as  the  array  argument  These  are  minimum  dimensions.  The 
arrays  can  be  dimensioned  larger  than  this  to  allow  for  program 
expansion. 

TABLE  2-1:  P.‘  J^ETERS  IN  THE  DOT  ARGUMENT  LIST 


PARAMETER 

DEFINITION 

INFO 

Information  parameter.  Before  calling  DOT 
the  first  time,  set  INFO-0.  When  control 
returns  from  DOT  to  the  calling  program, 

INFO  will  normally  have  a  value  of  0  or  1 . 

If  INFO-  0,  the  o^imization  is  complete  (or 
terminated  with  an  error  message).  If 

INFO-  1,  the  user  must  evaluate  the  objec¬ 
tive,  OBJ,  and  constraint  functions,  G(l), 
i-1  ,NCON,  and  call  DOT  again.  A  third 
possibility,  INFO-  2,  exists  also.  In  this 
case,  the  user  must  provide  gradient  in- 
fonnation.  This  is  an  advanced  feature 
and  is  described  in  Section  3.2. 

METHOD 

Optimization  method  to  be  used. 

METHOD  -  0  or  1  means  use  the  modified 
method  of  feasbie  directions. 

METHOD  -  2  means  use  the  sequential 
linear  programming  method.  If  the  problem 
in  unconstrained  (NCON  -  0)  the  BFGS  al¬ 
gorithm  will  be  used,  regardless  of  the 
value  of  the  parameter  METHOD. 
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DOT  WITH  APPLICATION  PROGRAMS 


IPRINT 

1 

Print  control  parameter. 

IPRINT  «  0  no  output. 

IPRINT  >  1  internal  parameters,  initial 
information  and  results. 

IPRINT  «  2  same  plus  objective  functton 
and  X-vector  at  each  iteration. 

IPRINT  B  3  same  plus  G-vector  and 

critical  constraint  numbers. 

IPRINT  -  4  same  plus  gradients. 

IPRINT  B  5  same  plus  search  direction. 

NOTE:  The  iPRM  Array  contains  additional 

print  options. 

NDV 

Number  of  decision/design  variables  con¬ 
tained  in  vector  X.  NDV  is  the  same  as  N 
in  the  mathematical  problem  statement 
given  in  Section  1.7. 

NCON 

Number  of  constraint  values  contained  in 
array  G.  NCON  is  the  same  as  M  in  the 
mathematical  problem  statement  given  in 
Section  1.7.  NCONbO  is  allowed. 

X(NDV) 

Vector  containing  the  design  variables.  On 
the  first  call  to  DOT.  this  is  the  user’s  best 
guess  of  the  design.  On  the  final  return 
from  DOT  (INFObO  is  returned),  the  vector 

X  contains  the  optimum  design. 

XL(NDV) 

Array  containing  lower  bounds  on  the 
design  variables,  X.  If  no  lower  bounds  are 
imposed  on  one  or  more  of  the  design  vari¬ 
ables,  the  corresportding  cornponent(s}  of 

XL  must  be  set  to  a  large  negative  number, 
say  -I.OE-t-IS.  Besureit’s-1.0E<fl5and 
not  -1.0E-15  (+15,  not  -15). 

2-6  2.2  PARAMETERS  IN  THE  CALLING  STATEMENT 
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DOT  WITH  APPLICATION  PROGRAMS 


XU(NDV)  Array  corrtaining  upper  bounds  on  the 
design  variables,  X.  If  no  upper  bounds 
are  imposed  on  one  or  more  of  the  design 
variables,  the  corresponding  component(s) 
of  XU  must  be  set  to  a  large  positive  num¬ 
ber,  say  1.0  E+15. 


Value  of  the  objective  function  correspond¬ 
ing  to  the  current  values  of  the  design  vari¬ 
ables  contained  in  X.  On  the  first  call  to 
DOT,  OBJ  need  not  be  defined.  DOT  will 
return  a  value  of  INFO-1  to  indicate  that 
the  user  must  evaluate  OBJ  and  call  DOT 
again.  Subsequently,  any  time  a  value  of 
INFO-1  is  returned  from  DOT,  the  objec¬ 
tive,  OBJ,  must  be  evaluated  for  the  cur¬ 
rent  design  and  DOT  must  be  called  again. 
OBJ  has  the  same  meaning  as  F(X)  in  the 
mathematical  problem  statement  given  in 
Section  1.7. 


MINMAX  Integer  parameter  specifying  whether  the 
minimum  {MINMAX-0,-1)  or  maximum 
(MINMAX-l)  of  the  objective  function  is  to 
be  found. 


G{NCON)  Array  containing  the  NCON  inequality  con¬ 
straint  values  corresponding  to  the  current 
design  contained  in  X.  On  the  first  call  to 
DOT,  the  constraint  values  need  not  be 
defined  On  return  from  DOT,  if  INFO-1, 
the  constraints  must  be  evaluated  for  the 
current  X  and  DOT  must  be  called  again. 

If  NCON-0,  array  G  must  be  dimensioned 
to  1  or  larger,  but  no  constraint  values 
need  to  be  provided. 


2.2  PARAMETERS  iN  THE  CALLING  STATEMENT 
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DOT  WITH  APPLICATION  PROGRAMS 


RPRM(20) 

Array  containing  the  real  (floating  point 
numbers)  control  parameters.  Initialize  the 
entire  array  to  0.0  to  use  all  default  values. 

If  you  use  other  values  than  the  defaults, 
set  the  corresponding  entries  to  the 
desired  values.  Section  3.1  describes  how 
to  change  the  value  of  these  parameters. 

IPRM(20) 

Array  containing  the  integer  control 
parameters.  As  with  the  RPRM  array,  set 
the  array  to  zero  to  use  the  default  values, 
or  set  the  proper  entries  to  the  desired 
values.  Section  3.1  describes  how  to 
change  the  value  of  these  parameters. 

WK(NRWK) 

User  provided  work  array  for  real  (ftoating 
point)  variables.  Array  WK  is  used  to  store 
internal  scalar  variables  and  arrays  used 
by  DOT.  if  the  user  has  not  provided 
enough  storage.  DOT  will  print  the  ap¬ 
propriate  message  and  terminate  the  op- 
tirr^tion. 

NRWK 

Dimensioned  size  of  work  anay  WK. 

NRWK  should  be  set  quite  large,  starting  at 
about  500  for  a  small  problem.  If  NRWK 
has  been  given  too  small  a  value,  an  error 
message  will  be  printed  and  the  optimiza¬ 
tion  will  be  terminated. 

IWK(NRIWK) 

User  provided  work  array  for  integer  (fixed 
point)  variables.  Array  IWK  is  used  to  store 
internal  scalar  variables  and  arrays  used 
by  DOT.  if  the  user  has  not  provided 
enough  storage,  DOT  will  print  the  ap¬ 
propriate  message  and  terminate  the  op¬ 
timization. 

2-8  2.2  PARAMETERS  IN  THE  CALLING  STATEMENT 
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DOT  WITH  APPLICATION  PROGRAMS 


NRIWK 

Dimensioned  size  of  work  array  IWK.  A 

good  estimate  is  300  for  a  small  problem. 

Increase  the  size  of  NRIWK  as  the  problem 

grows  larger.  If  NRIWK  is  too  small,  an 

error  message  will  be  printed  and  the  op- 

timization  will  be  terminated. 

Note:  The  minimum  required  values  of  NRWK  and  NRIWK  are 
defined  as  follows  (The  dimensions  may  be  larger  than  this): 

N1  =  NCON  -r  NDV 

N2  =  2*NDV 

N3  =  10*NDV 

N4  =  MIN  (Ni;^2) 

N5='l 


IF  NCON  =  0.  N5=0 
NCOLA=MAX(N3;^4) 

NGMAX  =  MIN  (NCON  J^COLA) 

NRB  =  MIN  (NlJ^COLA+1) 

IF  NCON  =  0,  NRB  =  1 

NRWK  =  NDV*(10+NCOLA)  +  5*NCON  +  NCOLA+  NRB**2  + 
MAX(NDV^COLA)+2*N5*NRB  +  40 

IF  METHOD  >  1,  NRWK  =  NRWK  +  NGMAX  + 

NDV*(3+NCOLA) 

NRIWK  =  NDV  +  NCON  +  NGMAX  +  71 
IF  METHOD  >  1.  NRIWK  =  NRIWK  +  NGMAX 


2.2  PARAMETERS  IN  THE  CALLING  STATEMENT  2-9 
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DOT  WITH  APPLICATION  PROGRAMS 


A  program  called  DTSTOR  is  provided  with  DOT.  If  you  com¬ 
pile  and  run  this  program  interactively,  the  minimum  required 
values  of  NRWK  and  NRIWK  are  calculated  for  you.  See  Ap¬ 
pendix  B  for  more  information  cn  this  option. 


2-10  2.2  PARAMETERS  IN  THE  CAUING  STATEMENT 
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APPENDIX  D 


VERIFICATION  AND  CASE  STUDY  OUTPUT 


VERIFICATION  #1 
OPTIMIZATION  SOLUTION 


A)  Problem  Parameters: 

Arch  Angle  :  0.003 

Arch  Radius:  1000000.000 
Arch  Height:  3.000 

E)  Derived  Constants: 

No  of  System  Nodal  Points... 
No  of  Degrees  of  Freedom. . . . 

Length  per  Element . 

Number  of  Iterations . 

C)  Structure  Loading: 

FX . 

FY . 

FM . 

FA . 


Youngs  Modulus:  30000000.0 
Yield  Strength:  52000.0 
No  of  Elements:  4 


5 

15 

11.2500 

0 


1000.0000 

0.0000 

0.0000 

0.0000 


Volume 

50.62481 

50.62481 

50.62481 

50.62481 

E)  Objective  Function: 

Total  structure  Volume:  202.499222 
Node  Stress 

1  19999.90 

2  14999.94 

3  9999.967 

4  4999.973 

5  9.7119728E-06 


D)  Elemental 

Dimensions 

and  Stress 

Distribution: 

Element 

Height 

Base 

Length 

1 

3.00000 

1.50000 

11.24996 

2 

3.00000 

1.50000 

11.24996 

3 

3.00000 

1.50000 

11.24996 

4 

3.00000 

1.50000 

11.24996 

F)  Boundary  Conditions: 
Node  X-Displ 

1  1 

5  0 


Y-Displ 

1 

0 


Slope 

1 

0 


G)  Solution  Vector: 
Node  X-Displ 

1  O.OOOOOOE+00 

2  0.257809E-01 

3  0.937488E-01 

4  0.189841E+00 

5  0.299996E+00 


Y-Displ 

O.OOOOOOE+00 

0.112328E-08 

0.409060E-08 

0.828730E-08 

0.130987E-07 


Slope 

O.OOOOOOE+00 

-0.437496E-02 

-0.749993E-02 

-0.937491E-02 

-0.999991E-02 
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VERIFICATION  #2 
OPTIMIZATION  SOLUTION 


A)  Problem  Parameters: 

Arch  Angle  :  0.003 

Arch  Radius:  1000000.000 
Arch  Height:  3.000 

B)  Derived  Constants: 

No  of  System  Nodal  Points... 
No  of  Degrees  of  Freedom. . . . 

Length  per  Element . 

Number  of  Iterations . 


Youngs  Modulus:  30000000.0 
Yield  Strength:  52000.0 
No  of  Elements :  4 


5 

15 

11.2500 

0 


C)  Structure  Loading: 


FX .  0.0000 

. .  1000.0000 

. .  0.0000 

FA .  0.0000 


D)  Elemental  Dimensions 


Element 

1 

2 

3 

4 


Height 

3.00000 

3.00000 

3.00000 

3.00000 


and  Stress  Distribution: 


Base 

1.50000 

1.50000 

1.50000 

1.50000 


Length 

11.24996 

11.24996 

11.24996 

11.24996 


Volume 

50.62481 

50.62481 

50.62481 

50.62481 


E)  Objective  Function: 

Total  structure  Volume:  202.499222 
Node  Stress 

1  222.2231 

2  222.2229 

3  222.2227 

4  222.2224 

5  222.2222 


F)  Boundary  Conditions: 
Node  X-Displ 

1  1 

5  0 


Y-Displ 

1 

0 


Slope 

1 

0 


G)  Solution  Vector: 
Node  X-Displ 

1  O.OOOOOOE+00 

2  0.112328E-08 

3  0.409060E-08 

4  0.828730E-08 

5  0.130987E-07 


Y-Displ 

O.OOOOOOE+00 

0.833330E-04 

0.166666E-03 

0.249999E-03 

0.333332L-03 


Slope 

O.OOOOOOE+00 

-0.191236E-09 

-0.327832E-09 

-0.409790E-09 

-0.437110E-09 
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VERIFICATION  #3 


OPTIMIZATION  SOLUTION 


A)  Problem  Parameters: 

Arch  Angle  :  0.003  Youngs  Modulus:  30000000.0 

Arch  Radius:  1000000.000  Yield  Strength:  52000.0 

Arch  Height:  3.000  No  of  Elements:  4 

B)  Derived  Constants: 

No  of  System  Nodal  Points... 

No  of  Degrees  of  Freedom. . . . 

Length  per  Element . 

Number  of  Iterations . 

C)  Structure  Loading: 


FX .  0.0000 

FY .  0.0000 

FM .  10000.0000 

FA .  0.0000 


D)  Elemental  Dimensions  and  Stress  Distribution: 

Element  Height  Base  Length 

1  3.00000  1.50000  11.24996 

2  3.00000  1.50000  11.24996 

3  3.00000  1.50000  11.24996 

4  3.00000  1.50000  11.24996 

E)  Objective  Function: 

Total  structure  Volume:  202.499222 
Node  Stress 

1  4444.438 

2  4444.438 

3  4444.441 

4  4444.441 

5  4444.441 

F)  Boundary  Conditions: 

Node  X-Displ  Y-Displ  Slope 

11  1  1 
5  0  0  0 

G)  Solution  Vector: 

Node  X-Displ  Y-Displ  Slope 

1  O.OOOOOOE+00  O.OOOOOOE+00  O.OOOOOOE+00 

2  -0.624994E-02  -0 . 273194E-09  O.lllllOE-02 

3  -0.249998E-01  -0 . 109277E-08  0.222221E-02 

4  -0.562495E-01  -0 . 245874E-08  0.333332E-02 

5  -0.999991E-01  -0 . 437110E-08  0.444442E-02 


Volume 

50.62481 

50.62481 

50.62481 

50.62481 
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VERIFICATION  #4 
OPTIMIZATION  SOLUTION 


A)  Problem  Parameters: 

Arch  Angle  :  90.000 

Arch  Radius:  45.000 

Arch  Height:  3.000 

B)  Derived  Constants: 

No  of  System  Nodal  Points... 
No  of  Degrees  of  Freedom. . . . 

Length  per  Element . 

Number  of  Iterations . 

C)  Structure  Loading: 

FX . 

FY . 

FM . 

FA . 


Youngs  Modulus:  30000000.0 
Yield  Strength:  52000.0 
No  of  Elements:  4 


5 

15 

17.5581 

0 


0.0000 

1000.0000 

0.0000 

0.0000 


Volume 
79.01159 
79.01159 
79.01159 
79.01159 

E)  Objective  Function: 

Total  structure  Volume:  316.046356 
Node  Stress 

1  20217.61 

2  18601.18 

3  14266.01 

4  7777.182 

5  43.39704 


D)  Elemental  Dimensions  and  Stress  Distribution: 
Element  Height  Base  Length 

1  3.00000  1.50000  17.55813 

2  3.00000  1.50000  17.55813 

3  3.00000  1.50000  17.55813 

4  3.00000  1.50000  17.55813 


F)  Boundary  Conditions: 
Node  X-Displ 

1  1 

5  0 


Y-Displ 

1 

0 


Slope 

1 

0 


G)  Solution  Vector: 
Node  X-Displ 

1  O.OOCOOOE+00 

2  -0.654617E-01 

3  -0.223502E+00 

4  -0.381543E+00 

5  -0.447005E+00 


Y-Displ 

O.OOOOOOE+00 

0.131512E-01 

0.118880E+00 

0.355536E+00 

0.684770E+00 


Slope 

O.OOOOOOE+00 

0.750656E-02 

0.138705E-01 

0.181227E-01 

0.196159E-01 
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OPTIMIZATION  #1 


OPTIMIZATION  SOLUTION 


A)  Problem  Parcuneters : 

Arch  Angle  :  0.002 

Arch  Radius:  1000000.000 
Arch  Height:  2.000 

B)  Derived  Constants: 

No  of  System  Nodal  Points... 
No  of  Degrees  of  Freedom. . . . 

Length  per  Element . 

Number  of  Iterations . 

C)  Structure  Loading: 

FX . 

FY . 

FM . 

FA . 


Youngs  Modulus:  30000000.0 
Yield  Strength:  52000.0 
No  of  Elements:  4 


5 

15 

6.0000 

0 


2000.0000 

0.0000 

0.0000 

0.0000 


Volume 
29.52806 
22.14510 
14.76169 
7.38456 

E)  Objective  Function: 

Total  structure  Volume:  73.819420 

Node  Stress 

1  52018.20 

2  52020.44 

3  52026.59 

4  52000.27 

5  9.4704286E-05 


D)  Elemental  Dimensions  and  Stress  Distribution: 


Element 

1 

2 

3 

4 


Height 

2.00000 

2.00000 

2.00000 

2.00000 


Base 

1.84550 

1.38407 

0.92261 

0.46154 


Length 

8.00000 

8.00000 

8.00000 

8.00000 


F)  Boundary  Conditions: 
Node  X-Displ 

1  1 

5  0 


Displ 

1 

0 


Slope 

1 

0 


G)  Solution  Vector: 
Node  X-Displ 

1  O.OOOOOOE+00 

2  0.508622E-01 

3  0.197286E+00 

4  0.433113E+00 

5  0.742914E+00 


Y-Displ 

O.OOOOOOE+00 

0.221694E-08 

0.860890E-08 

0.189046E-07 

0.324212E-07 


Slope 

O.OOOOOOE+00 

-0.121376E-01 

-0.236977E-01 

-0.341030E-01 

-0.410363E-01 
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OPTIMIZATION  #1 
OPTIMIZATION  SOLUTION 


A)  Problem  Parameters: 

Arch  Angle  :  0.002 

Arch  Radius:  1000000.000 
Arch  Height:  2.000 

B)  Derived  Constants: 

No  of  System  Nodal  Points... 
No  of  Degrees  of  Freedom. . . . 

Length  per  Element . 

Number  of  Iterations . 


Youngs  Modulus:  30000000.0 
Yield  Strength:  52000.0 
No  of  Elements:  8 


9 

27 

4.0000 

1 


2000.0000 

0.0000 

0.0000 

0.0000 


C)  structure  Loading: 

FX . . 

FY . 

FM . 

FA . 


D)  Elemental  Dimensions  and  Stress  Distribution: 


Element  Height 

Base 

Length 

Volume 

1 

2.00000 

1.84548 

4.00000 

14.76385 

2 

2.00000 

1.61572 

4.00000 

12.92572 

3 

2.00000 

1.38417 

4.00000 

11.07339 

4 

2.00000 

1.15351 

4.00000 

9.22805 

5 

2.00000 

0.92224 

4.00000 

7.37795 

6 

2.00000 

0.69289 

4.00000 

5.54310 

7 

2.00000 

0.46127 

4.00000 

3.69015 

8 

2.00000 

0.23665 

4.00000 

1.89324 

Objectiv®  Function; 

Total  structure  Volume:  66.495445 

Node 

Stress 

1 

52015.55 

2 

51986.30 

3 

52013.96 

4 

52013.29 

5 

52045.58 

6 

51955.51 

7 

52029.16 

8 

50706.71 

9 

0.2237021 

Boundary  Conditions: 

Node 

X-Displ 

Y-Displ 

Slope 

1 

1 

1 

1 

9 

0 

0 

0 

Solution  Vector: 

Node 

X-Displ 

Y-Displ 

Slope 

1 

O.OOOOOOE+00 

O.OOOOOOE+00 

O.OOOOOOE+00 

2 

0.132929E-01 

0.577893E-09  - 

-0.650196E-02 

3 

0.525036E-01 

0.228824E-08  - 

■0.129384E-01 

4 

0.117357E+00 

0.511886E-08  - 

■0.192957E-01 

5 

0.207485E+00 

0.905344E-08  - 

-0.255373E-01 

6 

0.322357E+00 

0.140683E-07  - 

•0.316093E-01 

7 

0.461109E+00 

0.201.^50E-07  - 

•0.373821E-01 

8 

0.622199E+00 

0.271538E-07  - 

■0.425851E-01 
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9  0.801554E+00 


0.349690E-07  -0 . 459655E-01 
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OPTIMIZATION  #1 
OPTIMIZATION  SOLUTION 


A)  Problem  Parameters: 

Arch  Angle  :  0.002 

Arch  Radius:  1000000.000 
Arch  Height:  2.000 

B)  Derived  Constants: 

No  of  System  Nodal  Points... 
No  of  Degrees  of  Freedom. . . . 

Length  per  Element . 

Number  of  Iterations . . 

C)  Structure  Loading: 

FX . 

FY . 

FM . 

FA . 


Youngs  Modulus: 

30000000.0 

Yield  Strength: 

52000.0 

No  of  Elements: 

12 

13 

39 

2.6667 

4 

2000.0000 

0.0000 

0.0000 

0.0000 

D)  Elemental  Dimensions  and  Stress  Distribution: 


Element 

Height 

Base 

Length 

Vo  1  lime 

1 

2.00000 

1.84829 

2.66667 

9.85756 

2 

2.00000 

1.69537 

2.66667 

9.04196 

3 

2.00000 

1.54275 

2.66667 

8.22798 

4 

2.00000 

1.38887 

2.66667 

7.40732 

5 

2.00000 

1.23566 

2.66667 

6.59016 

6 

2.00000 

1.08215 

2.66667 

5.77148 

7 

2.00000 

0.93133 

2.66667 

4.96712 

8 

2.00000 

0.77510 

2.66667 

4.13389 

9 

2.00000 

0.63162 

2.66667 

3.36863 

10 

2.00000 

0.46759 

2.66667 

2.49383 

11 

2.00000 

0.31955 

2.66667 

1.70426 

12 

2.00000 

0.20000 

2.66667 

1.06667 

E)  Objective  Function: 

Total  structure  Volume:  64.630867 

Node  Stress 

1  51925.73 

2  51892.89 

3  51843.36 

4  51329.66 

5  51784.52 

6  51739.70 

7  51531.16 

8  51598.32 

9  50656.20 

10  51320.27 

11  50069.02 

12  39998.66 

13  0.6707708 

F)  Boundary  Conditions: 

Node  X-Displ  Y-Displ  Slope 

11  1  1 
13  0  0  0 


G)  Solution  Vector: 

Node  X-Displ  Y-Displ  Slope 
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1  O.OOOOOOE+00 

2  0.598324E-02 

3  0.237427E-01 

4  0.532194E-01 

5  0.943462E-01 

6  0.147043E+00 

7  0.211205E+00 

8  0.286684E+00 

9  0.373299E+00 

10  0.470718E+00 

11  0.578546E+00 

12  0.696051E+00 

13  0.820672E+00 


O.OOOOOOE+00 

0.259434E-09 

0.103343E-08 

0.231938E-08 

0.411430E-08 

0.641459E-08 

0.921563E-08 

0.125107E-07 

0.162918E-07 

0.205440E-07 

0.252490E-07 

0.303731E-07 

0.358010E-07 


O.OOOOOOE+00 

-0.442333E-02 

-0.882640E-02 

-0.132043E-01 

-0.175555E-01 

-0.218709E-01 

-0.261415E-01 

-0.303404E-01 

-0.344682E-01 

-0.384082E-01 

-0.422098E-01 

-0.455477E-01 

-0.473255E-01 
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OPTIMIZATION  #2 


OPTIMIZATION  SOLUTION 


A)  Problem  Parameters: 

Arch  Angle  :  90.000  Youngs  Modulus:  30000000.0 

Arch  Radius:  32.000  Yield  Strength:  52000.0 

Arch  Height:  2.000  No  of  Elements:  12 

B)  Derived  Constants: 

No  of  System  Nodal  Points...  13 
No  of  Degrees  of  Freedom. ...  39 

Length  per  Element .  4.1858 

Number  of  Iterations .  1 

C)  Structure  Loadihg: 


FX .  0.0000 

FY .  -2000.0000 

FM .  0.0000 

FA .  0.0000 


D)  Elemental  Dimensions  and  Stress  Distribution: 


Element 

Height 

Base 

Length 

Volume 

1 

2.00000 

1.86405 

4.18580 

15.60506 

2 

2.00000 

1.84215 

4.18580 

15.42178 

3 

2.00000 

1.79464 

4.18580 

15.02403 

4 

2.00000 

1.71836 

4.18580 

14.38546 

5 

2.00000 

1.61213 

4.18580 

13.49611 

6 

2.00000 

1.47665 

4.18580 

12.36194 

7 

2.00000 

1.32535 

4.18580 

11.09533 

8 

2.00000 

1.16593 

4.18580 

9.76072 

9 

2.00000 

0.93161 

4.18580 

7.79905 

10 

2.00000 

0.71259 

4.18580 

5.96554 

11 

2.00000 

0.51683 

4.18580 

4.32666 

12 

2.00000 

0.36557 

4.18580 

3.06044 

E)  Objective  Function: 

Total  structure  Volume:  128.302124 
Node  Stress 

1  51975.23 

2  51991.84 

3  52001.67 

4  51952.61 

5  51914.72 

6  51926.68 

7  51573.75 

8  50482.76 

9  51888.91 

10  51924.89 

11  48451.66 

12  34653.70 

13  178.9231 

F)  Boundary  Conditions: 

Node  X-Displ  Y-Displ  Slope 

11  1  1 
13  0  0  0 

G)  Solution  Vector: 

Node  X-Displ  Y-Displ  Slope 
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1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 


O.OOOOOOE+00 

0.149422E-01 

0.589253E-01 

0.129477E+00 

0.222644E+00 

0.333172E+C0 

0.454693E+00 

0.579860E+00 

0.700496E+00 

0.808353E+00 

0.895118E+00 

0.952352E+00 

0.972583E+00 


O.OOOOOOE+00 
-0.105421E-02 
-0.987875E-02 
-0.339057E-01 
-0.799319E-01 
-0.153870E+00 
-0.260536E+00 
-0.403368E+00 
-0.58^032E+00 
-0.802894E+00 
-0. 105869E+01 
-0.134670E+01 
-0.165575E+01 


O.OOOOOOE+00 

-0.714707E-02 

-0.142563E-01 

-0.213028E-01 

-0.282739E-01 

-0.351633E-01 

-0.419651E-01 

-0.486117E-01 

-0.549785E-01 

-0.613227E-01 

-0.673512E-01 

-0.723963E-01 

-0.747875E-01 
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OPTIMIZATION  #3 


OPTIMIZATION  SOLUTION 


A)  Problem  Parameters: 

Arch  Angle  :  90.000  Youngs  Modulus:  30000000.0 

Arch  Radius:  32.000  Yield  Strength:  52000.0 

Arch  Height:  2.000  No  of  Elements:  12 

B)  Derived  Constants: 

No  of  System  Nodal  Points...  13 
No  of  Degrees  of  Freedom. ...  39 

Length  per  Element .  4.1858 

Number  of  Iterations .  1 

C)  Structure  Loading: 


FX .  2000.0000 

FY .  0.0000 

FM .  0.0000 

FA .  0.0000 


D)  Elemental  Dimensions  and  Stress  Distribution: 


Element 

Height 

Base 

Length 

Volume 

1 

2.00000 

1.84614 

4.18580 

15.45518 

2 

2.00000 

1.68690 

4.18580 

14.12203 

3 

2.00000 

1.43532 

4.18580 

12.01594 

4 

2.00000 

1.19577 

4.18580 

10.01052 

5 

2.00000 

0.97029 

4.18580 

8.12291 

6 

2.00000 

0.75740 

4.18580 

6.34069 

7 

2.00000 

0.56761 

4.18580 

4.75181 

8 

2.00000 

0.40093 

4.18580 

3.35641 

9 

2.00000 

0.26110 

4.18580 

2.18585 

10 

2.00000 

0.20000 

4.18580 

1.67432 

11 

2.00000 

0.20000 

4.18580 

1.67432 

12 

2.00000 

0.20000 

4.18580 

1.67432 

E)  Objective  Function: 

Total  structure  Volume:  81.384300 

Node  Stress 

1  52001.90 

2  51862.35 

3  51955.02 

4  51942.95 

5  51853.82 

6  51974.19 

7  51925.02 

8  51871.48 

9  51649.30 

10  41477.37 

11  21322.36 

12  9098.637 

13  4994.856 

F)  Boundary  Conditions: 

Node  X-Displ  Y-Displ  Slope 

11  1  1 
13  0  0  0 

G)  Solution  Vector: 

Node  X-Displ  Y-Displ  Slope 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 


O.OOOOOOE+00 

0.144840E-01 

0.557774E-01 

0.120917E+00 

0.206307E+00 

0.307129E+00 

0.417585E+00 

0.531063E+00 

0.640343E+00 

0.737912E+00 

0.813769E+00 

0.860985E+00 

0.877143E+00 


O.OOOOOOE+00 

-0.944371E-03 

-0.914169E-02 

-0.312207E-01 

-0.732726E-01 

-0.140544E+00 

-0.237250E+00 

-0.366366E+00 

-0.529394E+00 

-0.726163E+00 

-0.948447E+00 

-0.118418E+01 

-0.142556E+01 


O.OOOOOOE+00 

-0.677752E-02 

-0.131680E-01 

-0.195029E-01 

-0.257570E-01 

-0.319051E-01 

-0.379512E-01 

-0.438425E-01 

-0.495298E-01 

-0.549162E-01 

-0.570368E-01 

-0.577029E-01 

-0.578414E-01 


-85- 


OPTIMIZATION  #4 
OPTIMIZATION  SOLUTION 


A)  Problem  Parameters: 

Arch  Angle  :  90.000 

Arch  Radius:  32.000 

Arch  Height:  2.000 

B)  Derived  Constants: 

No  of  System  Nodal  Points . . . 
No  of  Degrees  of  Freedom. . . . 

Length  per  Element . 

Number  of  Iterations . 

C)  Structure  Loading: 

FX . 

FY . 

FM . 

FA . 


Youngs  Modulus: 

30000000.0 

Yield  Strength: 

52000.0 

No  of  Elements: 

12 

13 

39 

4.1858 

1 

0.0000 

0.0000 

0.0000 

-100.0000 

D)  Elemental  Dimensions  and  Stress  Distribution: 


Element  Height 

Base 

Length 

Volume 

1 

2.00000 

2.08018 

4.18580 

17.41440 

2 

2.00000 

1.93027 

4.18580 

16.15947 

3 

2.00000 

1.76172 

4.18580 

14.74842 

4 

2.00000 

1.55158 

4.18580 

12.98917 

5 

2.00000 

1.31450 

4.18580 

11.00444 

6 

2.00000 

1.06716 

4.18580 

8.93384 

7 

2.00000 

0.82127 

4.18580 

6.87536 

8 

2.00000 

0.59123 

4.18580 

4.94951 

9 

2.00000 

0.39007 

4.18580 

3.26550 

10 

2.00000 

0.22180 

4.18580 

1.85686 

11 

2.00000 

0.20000 

4.18580 

1.67432 

12 

2.00000 

0.20000 

4.18580 

1.67432 

Objective  Function: 

Total 

structure  Volume 

:  101.545610 

Node 

Stress 

1 

52103.21 

2 

52039.99 

3 

52033.59 

4 

52017.27 

5 

52015.93 

6 

51958.49 

7 

51950.00 

8 

51920.60 

9 

51610.81 

10 

51626.45 

11 

25188.80 

12 

5614.041 

13 

1096.744 

Boundary  Conditions : 

Node 

X-Displ  Y 

-Displ 

Slope 

1 

1 

1 

1 

13 

0 

0 

0 

G)  Solution  Vector: 

Node  X-Displ  Y-Displ  Slope 
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1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 


O.OOOOOOE+00 

0.145714E-01 

0.573794E-01 

0.125999E+00 

0.216392E+00 

0.323304E+00 

0.440430E+00 

0.560627E+00 

0.676171E+00 

0.778995E+00 

0.861090E+00 

0.913642E+00 

0.931478E+00 


O.OOOOOOE+00 

-0.111348E-02 

-0,978175E-02 

-0.332242E-01 

-0.779492E-01 

-0.149537E+00 

-0.252410E+00 

-0.389635E+00 

-0.562738E+00 

-0.771440E+00 

-0.101350E+01 

-0.127777E+01 

-0.154990E+01 


O.OOOOOOE+00 

-0.69174,tE-02 

-0.138459E-01 

-0.206572E-01 

-0.273459E-01 

-0.339032E-01 

-0.403032E-01 

-0.465215E-01 

-0.525061E-01 

-0.581423E-01 

-0.633131E-01 

-0.650579E-01 

-0.651839E-01 
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OPTIMIZATION  #5 
OPTIMIZATION  SOLUTION 


A)  Problem  Parameters: 

Arch  Angle  :  90,000 

Arch  Radius:  32.000 

Arch  Height:  2.000 

B)  Derived’ Constants: 

No  of  System  Nodal  Points... 
No  of  Degrees  of  Freedom. . . . 

Length  per  Element . 

Number  of  Iterations . 

C)  Structure  Loading: 

FX . 

FY . 

FM . 

FA . 


Youngs  Modulus:  30000000.0 
Yield  Strength:  52000.0 
No  of  Elements:  12 


13 

39 

4.1858 

1 


0.0000 

-8000.0000 

0.0000 

0.0000 


D) 


E) 


Elemental  Dimensions  and  Stress  Distribution: 


lement 

Height 

Base 

Length 

Volume 

1 

2.00000 

0,51296 

4.18580 

4.29430 

2 

2.00000 

0.90694 

4.18580 

7.59253 

3 

2.00000 

1.11658 

4.18580 

9,34755 

4 

2.00C00 

1.47904 

4.18580 

12.38193 

5 

2.00000 

1.19599 

4.18580 

10.01239 

6 

2.00000 

1.12703 

4.18580 

9.43503 

7 

2.00000 

1.53118 

4.18580 

12.81847 

8 

2,00000 

0.57650 

4.18580 

4.82624 

9 

2.00000 

0.59349 

4.18580 

4.96845 

10 

2.00000 

1.34594 

4.18580 

11.26764 

11 

2.00000 

2.20595 

4.18580 

18,46734 

12 

2.00000 

3.12364 

4.18580 

26.14985 

Objective  Function: 

Total 

structure  Volume: 

131.561722 

Node  Stress 


1 

8078.639 

2 

52000.93 

3 

50895.98 

4 

52005.46 

5 

52048.33 

6 

52010.10 

7 

42175,72 

8 

49131.36 

9 

4933.140 

10 

52013.68 

11 

52043.57 

12 

52041.75 

13 

52027.50 

F)  Boundary  Conditions: 

Node  X-Displ  Y-Displ 

1  1  1 

13  1  0 


Slope 

0 

1 


G)  Solution  Vector: 

Node  X-Displ  Y-Displ  Slope 
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1  O.OOOOOOE+00 

2  -0.823069E-01 

3  -0.143411E+00 

4  -0.178618E+00 

5  -0.189389E+00 

6  -0.177714E+00 

7  -0.146174E+00 

8  -0.105000E+00 

9  -0.619691E-01 
10  -0.273824E-01 

I  -0.803875E-02 

12  -0.764124E-03 

13  O.OOOOOOE+00 


O.OOOOOOE+00 

0.426506E-02 

0.157328E-01 

0.270851E-01 

0.319111E-01 

0.234612E-01 

-0.494672E-02 

-0.525035E-01 

-0.118718E+00 

-0.190906E+00 

-0.249018E+00 

-0.286585E+00 

-0.300012E+00 


0.208648E-01 

0-173341E-01 

0.118840E-01 

0.554722E-02 

0.244944E-04 

-0.683719E-02 

-0.132138E-01 

-0.165616E-01 

-0.201392E-01 

-0.168476E-01 

-0.117570E-01 

-0.606703E-02 

O.OOOOOOE-t-00 


OPTIMIZATION  #6 


OPTIMIZATION  SOLUTION 


A)  Problem  Parameters: 

Arch  Angle  :  90.000 

Arch  Radius:  32.000 

Arch  Height:  2.000 

B)  Derived  Constants: 

No  of  System  Nodal  Points... 
No  of  Degrees  of  Freedom. . . . 

Length  per  Element . 

Number  of  Iterations . 

C)  Structure  Loading: 

FX . 

FY . 

FM . 

FA . 


Youngs  Modulus:  30000000.0 
Yield  Strength:  52000.0 
No  of  Elements:  12 


13 

39 

4.1858 

1 


0.0000 

-8000.0000 

0.0000 

0.0000 


D)  Elemental  Dimensions  and  Stress  Distribution: 


Element  Height 

Base 

Length 

Volume 

1 

2.00000 

1.72707 

4.18580 

14.45838 

2 

2.00000 

0.94577 

4.18580 

7.91762 

3 

2.00000 

0.27194 

4.18580 

2.27661 

4 

2.00000 

0.65040 

4.18580 

5.44488 

5 

2.00000 

1.09075 

4.18580 

9.13130 

€ 

2.00000 

0.78986 

4.18580 

6.61236 

7 

2.00000 

1.06223 

4.18580 

8.89255 

8 

2.00000 

0.75388 

4.18580 

6.31116 

9 

2.00000 

0.25849 

4.18580 

2.16400 

10 

2.00000 

0.93275 

4.18580 

7.80861 

11 

2.00000 

1.75263 

4.18580 

14.67236 

12 

2.00000 

2.65143 

4.18580 

22.19675 

Objective  Function: 

Total 

structure  Volume 

:  107.886574 

Node 

Stress 

1 

51991.36 

2 

49986.21 

3 

46231.32 

4 

51991.33 

5 

49404.85 

6 

51991.04 

7 

51841.79 

8 

42080.21 

9 

51997.34 

10 

51279.07 

11 

52434.15 

12 

51908.94 

13 

51935.26 

Boundary  Conditions: 

Node 

X-Displ  Y 

-Displ 

Slope 

1 

1 

1 

1 

13 

1 

0 

1 

G)  Solution  Vector: 

Node  X-Displ  Y-Displ  Slope 


1  O.OOOOOOE+00 

2  -0.121854E-01 

3  -0.442765E-01 

4  -0.861533E-01 

5  -0.112986E+00 

6  -0.121058E+00 

7  -0.110454E+00 

8  -0.850803E-01 

9  -0.537900E-01 

10  -0.252856E-01 

11  -0.744947E-02 

12  -0.654979E-03 

13  O.OOOOOOE+00 


O.OOOOOOE+00 
0.457452E-03 
0.615066E-02 
0.177060E-01 
0.297132E-01 
0.342979E-01  - 
0.237526E-01  - 
-0.622361E-02  - 
-0.547404E-01  - 
-0.118485E+00  - 
-0.173190E+00  - 
-0.209063E+00  - 
-0.222069E+00 


O.OOOOOOE+00 

0.527930E-02 

0.947556E-02 

0.907372E-02 

0.435047E-02 

0.109393E-03 

0.701848E-02 

0.115516E-01 

0.155000E-01 

0.155501E-01 

0.110976E-01 

0.580765E-02 

O.OOOOOOE+00 


APPENDIX  E 

ARCH.OPTJOR  SOURCE  CODE 

PROGRAM  ARCH_OPTIMIZATION 

★★***★★**★**★**★***★★*****★*******★**★♦♦*★*********★****★*★★******★*★* 
★  * 

*  ARCH  OPTIMIZATION  ANALYSIS  CODE  * 

*  * 
********************************************************************** 
* 

*  ALPHA _ TRANSFORMATION  ANGLE  OF  ELEMENT  (ANGLE  TO  X-AXIS) 

*  ANGLE. .TOTAL  ANGLE  OF  ARCH  (IN  DEGREES) 

*  BASE . DOT  ARRAY  CONTAINING  THE  ELEMENTAL  BASE  DIMENSIONS 

*  BASEL.... DOT  ARRAY  CONTAINING  THE  ELEMENTAL  BASE  DIMENSIONS  LOWER 

*  SIDE  CONSTRAINT 

*  BASEU....DOT  ARRAY  CONTAINING  THE  ELEMENTAL  BASE  DIMENSIONS  UPPER 

*  SIDE  CONSTRAINT 

*  BETA _ TRANSFORMATION  ANGLE  OF  ELEMENT  (ANGLE  TO  Y-AXIS) 

*  B_1 . BOUNDARY  TERMS  APPLIED  AT  END  ”1" 

*  B_2 . BOUNDARY  TERMS  APPLIED  AT  END  “2" 

*  Cl C5 ...  CONSTANTS  RELATED  TO  ELEMENT  STIFFNESS  COEFFICIENTS 

*  CLAN . CONCENTRATED  LOAD  APPLICATION  NODE  (THE  NODE  FX,FY,FM  ARE 

*  APPLIED 

*  COUNT. .. .COUNTS  THE  NUMBER  OF  ITERATIONS  COMPLETED 

*  DOF . DEGREE  OF  FREEDOMS  (UNKNOWN  DISPLACEMENTS  &  SLOPES) 

*  DSN . DESIGN  VARIABLE  FOR  EACH  ELEMENT 

*  DVIBG. .. .DESIGN  VARIABLE  #1  (BASE  DIMENSION)  INITIAL  ESTIMATE 

*  DVILO. .. .DESIGN  VARIABLE  #1  (BASE  DIMENSION)  LOWER  SIDE  CONSTRAINT 

*  DVIUP. .. .DESIGN  VARIABLE  #1  (BASE  DIMENSION)  UPPER  SIDE  CONSTRAINT 

*  EK . 6X6  ELEMENT  STIFFNESS  MATRIX  IN  LOCAL  X, Y  COORDINATES 

*  EKPR . 6X6  ELEMENT  STIFFNESS  MATRIX  IN  ELEMENT  LOCAL  COORDINATES 

*  ELEN . LENGTH  OF  ELEMENT 

*  F . FORCE  VECTOR  OF  SYSTEM 

*  FA . CONSTANT  DISTRIBUTED  LOAD  IN  X  DIRECTION  FROM  END  TO  END 

*  FM . CONCENTRATED  MOMENT  AT  FREE  END 

*  FX . CONCENTRATED  LOAD  IN  X  DIRECTION  AT  FREE  END 

*  FY . CONCENTRATED  LOAD  IN  Y  DIRECTION  AT  FREE  END 

*  G . THE  ARRAY  OF  CONSTRAINT  FUNCTIONS 

*  GAMMA.... 6X6  ELEMENT  TRANSFORMATION  MATRIX 

*  GK . (NDOF)X(NDOF)  GLOBAL  STIFFNESS  MATRIX 

*  H . DEPTH  OF  ARCH  SECTION 

*  INDSN. .. .INITIAL  (UNIFORM)  DESIGN  DIMENSION 

*  INFO . DOT  PARAMETER  USED  TO  SIGNAL  THAT  THE  OPT  IS  COMPLETE 

*  IPRINT...DOT  PARAMETER  USED  SELECT  THE  DATA  OUTPUT  FORMAT 

*  IPRM . DOT  SELECTABLE  INTEGER  PARAMETERS 

*  ITERATE.. THE  NUMBER  OF  TIMES  DOT  IS  TO  BE  RELOADED  WITH  THE 

*  PRECEEDING  DATA 

*  IWK . DOT  INTERNAL  WORK  SPACE  ARRAY 

*  METHOD... DOT  PARAMETER  USED  TO  DEFINE  THE  OPTIMIZATION  METHOD 

*  MINMAX...DOT  PARAMETER  USED  TO  MINIMIZE/MAXIMIZE  THE  PROBLEM 

*  NCON . NUMBER  OF  DESIGN  CONSTRAINTS 

*  NDOF . NUMBER  OF  DEGREES  OF  FREEDOM 

*  NEL . NUMBER  OF  ELEMENTS 

*  NRIWK....DOT  INTERNAL  WORK  SPACE  ARRAY  DIMENSION 

*  NRWK . DOT  INTERNAL  WORK  SPACE  ARRAY  DIMENSION 


-92 


onooono  *»  o  no  on  on  on  no  on 


*  NSNP . NUMBER  OF  SYSTEM  NODAL  POINTS 

*  OBJ . THE  OBJECTIVE  FUNCTION  OF  THE  OPTIMIZATION 

*  OPTDCS. . .OPTIMIZATION  DECISION  TO  OPTIMIZE  THE  PROBLEM  OR  NOT 

*  PI. . .P4. .PARAMETER  DIMENSION  CORRESPONDING  TO  THE  NEL,  NSNP,  NCON, 

*  AND  NDOF,  RESPECTIVELY 

*  PHI . SUBTENDED  ELELENT  ANGLE  (ALSO,  PHIANG  IN  DEGREES) 


*  PRCSN....THE  PRECISION  DESIRED  TO  SOLVE  THE  FEM  SYSTEM  OF  EQUATIONS 

*  RADIUS. . .ARCH  RADIUS 

*  RPRM . DOT  SELECTABLE  REAL  PARAMETERS 

*  SIGMA_B..THE  ELEMENTAL  NORMAL  STRESS  DUE  TO  BENDING 

*  SIGMA_N..THE  ELEMENTAL  NORMAL  STRESS  DUE  TO  AXIAL  FORCES 

*  SIGMA_T..THE  MAXIMUM  TOTAL  STRESS  IN  EACH  ELEMENT 

*  U . THE  "DISPLACEMENT"  VECTOR  OF  THE  SYSTEM  OF  LINEAR  EQUATIONS 

*  WK . DOT  INTERNAL  WORK  AREA 

*  X . GLOBAL  HORIZONTAL  COORDINATE 

*  Y . GLOBAL  VERTICAL  COORDINATE 

*  YIELD. .. .YIELD  STRENGTH  OF  THE  ARCH  MATERIAL 

*  YOUNG. .. .YOUNG'S  MODULUS  OF  THE  ARCH  MATERIAL 

C  ....declare  the  variables . 

INCLUDE  ' ARCH_COM.FOR' 

....read  the  input  parameters . 

OPEN(8,  FILE-'ARCH_IN.DAT' ,  STATUS-' OLD' ) 

READ(8, *)  ANGLE, RADIUS, YOUNG, YIELD, NEL,  METHOD,  IPRINT,  DVIBG, 

&  DVILO, DVIUP, H, CLAN, FX, FY, FM, FA, OPTDCS, ITERATE, PRCSN, 

&  BX1,BY1,BM1,BX2,BY2,BM2, LABEL 

....define  constants . 

NCON  -  3*NEL 
NSNP  -  NEL  +  I 
NDOF  -  3*NSNP 

....determine  the  system  nodal  coord  and  element  orientation.. 
CALL  GEOMETRY  (NEL, NSNP, ANGLE, RADIUS, X, Y, ALPHA, BETA, ELEN) 

....define  the  size  of  the  worJc  arrays  for  DOT . 

NRWK  -  38800 
NRIWK  -  400 

....Optimize  the  problem . 

CALL  OPTIMIZATION_TOOL 

....con^ile  and  format  the  output . 

CALL  ARCH_OUTPUT 

END 

********************************************************************** 


SUBROUTINE  GEOMETRY  (NEL, NSNP, ANGLE, RADIUS, X, Y, ALPHA, BETA, ELEN) 


I  This  routine  is  used  by  main  ARCH_OPTIMIZATION  to  generate  I 
I  the  X-,  y-coordinates  of  each  system  node,  to  determine  I 
I  the  orientation  of  each  element,  and  to  calculate  the  I 
I  length  of  each  element.  I 


....declare  the  variables . 

INTEGER  NEL, NSNP, PI, P2 
PARAMETER  (Pl-32, P2-33) 

REAL  ANGLE, RADIUS, ELEN, X(P2) ,Y(P2) ,ALPHA(P1)  BETA(Pl), 
£  PI,PHI,ANG,YNUM,XDEN 


-93- 


oo  oo  oo  no  on  no  on  onoooo  »*  o  on  o  o  no 


PARAMETER (PI  -  3.141593) 

....determine  the  geometric  constants... 
PHI  -  (ANGLE/NEL)* (PI/180.0) 

X(l)  -  0.0 
T(1)  -  0.0 

DO  100  i-2,  NSNP 

ANG  -  (i-1.0)*PHI 
X(i)  -  RADIUS  *  (1.0  -  COS (ANG)) 

Y(i)  -  RADIUS  *  SIN (ANG) 

YNUM  -  (Y(i)  -  Y(i-l)) 

XDEN  -  (X(i)  -  X(i-l)) 

ALPHA(i-l)  -  AT AN2 (YNUM, XDEN) 
BETA(i-l)  -  (PI/2.0)  -  ALPHA(i-l) 

100  CONTINUE 

....determine  the  length  of  each  element 
ELEN  -  SORT (X (2) **2.0  +  Y(2)**2.0) 

RETURN 

END 


SUBROUTINE  OPTIMIZATION  TOOL 


This  subroutine  directs  the  program  flow  optimization  decision 
i.e.,  optimize  the  problem  or  not.  It  also  serves  to  set  up  £ 
execute  the  DOT  optimization  software. 

....declare  the  variables . 

INCLUDE  'ARCH  COM. FOR' 

INTEGER  i 

....zero  out  the  RPRM  and  IPRM  arrays . 

DO  100  i-1,20 
RPRM(i)  -  0.0 
IPRM(i)  -  0 
100  CONTINUE 

....initialize  COUNT . 

COUNT  -  0 

....refine  the  constraint  tolerence . 

RPRM(2)  -  0.0001 
RPRM(3)  -  0.0001 

....turn  off  DOT'S  auto  scaling . 

IPRM(2)  -  -1 

....increase  DOT'S  default  number  of  iterations . 

IPRM(3)  -  1000 
IPRM(8)  -  1000 

....increase  DOT'S  number  of  consecutive  convergence  criteria. 
IPRM(4)  -  3 
IPRM(9)  -  3 

....define  MINMAX— 1  to  minimize  the  objective  function . 

MINMAX  -  -1 
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....initialize  the  design  variable  limits  and  best  guess . 

DO  200  i-l,NEL 

BASE(i)  -  DVIBG 
BASEL <i)  -  DVILO 
BASED (i)  -  DVIUP 

200  CONTINUE 

. . . .make  optimization  decision . 

IF  (OPTDCS  .NE.  1)  THEN 
CALL  EVAL 
RETURN 
ENDIF 

....ready  to  optimize . 

INFO  -  0 

300  CALL  DOT  (INFO, METHOD,  IPRINT, NEL, NCON, BASE, BASEL, BASED, OBJ, 

&  MINMAX, G, RPRM, IPRM, WK, NRWK, IWK, NRIWK) 

....evaluate  the  objective  function  and  constraints . 

IF  (INFO  .GT.  0)  THEN 
CALL  EVAL 
GOTO  300 
ENDIF 

....refine  the  solution  vector  by  reoptimizing . 

IF  (COUNT  .LT.  ITERATE)  THEN 
INFO  -  0 
COUNT  -  COUNT +1 
GOTO  300 
ENDIF 

RETURN 

END 

******************************************************************** 

* 

SUBROUTINE  EVAL 


This  subroutine  is  used  to  evaluate  the  Objective  function, 
constraint  functions,  and  side  constraints  of  the  optimization 
problem. 


....declare  the  variables . 

INCLUDE  'ARCH_COM.FOR' 

INTEGER  i, j 

....calculate  the  objective  function . 

OBJ  -  0.0 

DO  100  i-l,NEL 

OBJ  -  OBJ  +  BASE(i)*H*ELEN 
100  CONTINUE 

....initialize  the  design  constraint  vector 
DO  200  i-l,NCON 
G(i)  -  0.0 
200  CONTINUE 

....determine  the  design  constraints . 

CALL  ARCH  STRESS 
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C 

DO  210  i“l,NEL 

IF  (SIGMA_T(i)  .GE.  SIGMA_T (i+1) )  THEN 
SIGMA  -  SIGMA_T(i) 

ELSE 

SIGMA  -  SIGMA_T(i+l) 

ENDIF 

C 

G(i)  -  (SIGMA/YIELD)  -  1.0 
210  CONTINUE 
C 

DO  220  i»l,:’EL 
j-=i+NEL 

G(j)  -  (BASE(i) / (3.0*H) )  -  1.0 
220  CONTINUE 
C 

DO  230  i-l,NEL 

*i+(2*NEL) 

G(j)  -  H/ (l0.0*BASE(i) )  -  1.0 
230  CONTINUE 
C 

RETURN 

END 

********************************************************************** 


* 


SUBROUTINE  ARCH  STRESS 


It 


This  subroutine  is  used  to  perform  the  Finite  Element  analysis 
of  the  stresses  developed  in  an  arch  or  beam  for  a  given  load¬ 
ing. 


....declare  the  variables . 

INCLUDE  'ARCH_COM.FOR' 

INTEGER  I PVT (99) 

REAL  GK(P4,P4) ,F(P4) 

REAL*8  BK<P4,P4) ,BF(P4) ,BU(P4) ,FAC(9801) ,W0RK(99) 


....form  the  element  and  system  matrices . 

CALL  FORM  (NEL, NDOF, ALPHA, BETA, H, ELEN, YOUNG, BASE, GK) 

....form  the  Force  vector,  F . 

CALL  FORCE_VECTOR  (NEL, NDOF, ELEN, ALPHA, BETA, FA, F) 

....set  the  boundary  conditiona  and  loads . 

CALL  BNDARY  (NDOF,  GK,  CLAN,  FX,  FY,  FM,  F,  BXl ,  BYl ,  BMl,  BX2,  BY2,  BM2) 

....solve  the  system  oi  equations . 

IF  (PRCSN  .EQ.  2)  THEN 

. . . .change  GK  and  F  arrays  to  double  precision. ...... 

CALL  UPSCALE  (NDOF, GK, F, BK,  BF) 

....solve  the  system  of  equations . 

CALL  DL2ARG  (NDOF, BK, P4 , BF, 1, BU, FAC, IPVT, WORK) 

....change  BU  array  to  single  presicion . 

CALL  DOWNSCALE  (NDOF,BU,U) 

ELSE 

....solve  the  system  of  equations . 

CALL  L2ARG  (NDOF, GK, P4 , F, 1 , U, FAC, IPVT, WORK) 

ENDIF 

....determine  the  stress  distribution . 
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CALL  STRESS (X, Y, ALPHA, BETA, U, NEL, ELEN, YOUNG, H, SIGMA_T) 

C 

RETURN 

END 

***************************************'^****************************** 
»  * 
SUBROUTINE  FORM  (NEL, NDOF, ALPHA, BETA, H, ELEN, YOUNG, BASE, GK) 


This  subroutine  is  used  to  construct  the  global  stiffness  mat¬ 
rix  for  the  arch  problem. 


....declare  the  variables . 

INTEGER  NEL, NDOF, NPOW, lEL, I , J, K, II , JJ, KK, I II , JJJ, PI , P2 , P4 

PARAMETER (Pl-32, P2-33, P4-99) 

REAL  ELEN, H, BASE (PI) , ALPHA (PI ), BETA (PI ) , YOUNG, 

&  C1,C2,C3,C4,C5,CA,CB,EKPR(6, 6) , GAM (6, 6) ,EK(P1, 6, 6) , 

£  GAMMA(P1,€,6) ,GK(P4,P4),EKGA(6,6) ,GAEKGA(6,6) , 

&  ALPHAI , BETAI 

. . . .define  the  constants  Cx . 

NPOW  -  1 

Cl  -  YOUNG* H/ ELEN 

C2  -  (H/ELEN) **2.0 

C3  -  (H**2.0) / (2.0*ELEN) 

C4  -  (H**2.0)/3.0 
C5  -  C4/2.0 

....initialize  the  EKPR  and  GAM  arrays . 

DO  100  II  -  1,6 
DO  90  JJ-  1,6 

EKPR (II, JJ)  -  O.U 
GAM(II,JJ)  -  0.0 
90  CONTINUE 

100  CONTINUE 

....initialize  the  EK  and  GAMMA  arrays . 

DO  130  lEL  -  1,NEL 
DO  120  1-1,6 
DO  110  J  -  1,6 

EK(IEL,1,J)  -  0.0 
GANMAdEL,  I,  J)  -  0.0 
110  CONTINUE 

120  CONTINUE 

130  CONTINUE 


....determine  the  EKPR 

EKPR (1, 1) 

» 

Cl 

EKPRd,  4) 

- 

-Cl 

EKPR(2,2) 

- 

C1*C2 

EKPR(2, 3) 

■ 

C1*C3 

EKPR(2, 5) 

m 

-C1*C2 

EKPR(2,  6) 

m 

C1*C3 

EKPR(3,2) 

- 

C1*C3 

EKPR(3,3) 

- 

C1*C4 

EKPR(3, 5) 

-C1*C3 

EKPR(3, 6) 

- 

C1*C5 

EKPR (4, 1) 

- 

-Cl 

EKPR(4,4) 

- 

Cl 

EKPR(5,2) 

- 

-C1*C2 
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EKPR(5, 3) 
EKPR(5, 5) 
EKPR(5,  6) 
EKPR(6, 2) 
EKPR(6,3) 
EKPR(6,5) 
EKPR(6, 6) 


-  -C1*C3 

-  C1*C2 

-  -C1*C3 

-  C1*C3 

-  C1*C5 

-  -C1*C3 

-  C1*C4 


....initialize  the  GK  array 
DO  150  I  -  1,  NDOF 
DO  140  J  -  1,  NDOF 
GK(1,J)  -  0.0 
140  CONTINUE 

150  CONTINXre: 


. . . .determine  the  GAMMA  matrix 
DO  170  lEL  -  1,NEL 

ALPHAI  -  ALPHA (lEL) 

BETAI  -  BETA(IEL) 

CA  -  COS (ALPHAI) 

CB  -  COS (BETAI) 

CaVMMAdEL,  1, 1)  -  CA 
GAMMA (lEL, 1,2)  -  CB 
GAMMA (lEL, 2,1)  -  -CB 
GAMMA(IEL,2,2)  -  CA 
GAMMA(IEL,3,3)  -  1.0 

GAMMA(IEL,4,4)  -  CA 
GAMMA(IEL,4,5)  -  CB 
GAMMAdEL,  5,  4)  -  -CB 
GAMMA(IEL,5,5)  -  CA 
GAMMAdEL,  6,  6)  -  1.0 

170  CONTINUE 


....  initialize  the  EKGA  and  GAEKGA  arrays 
DO  270  lEL  -  1,  NEL 
DO  190  III  -  1,6 
DO  180  JJJ  -  1,6 

EKGA (III, JJJ)  -  0.0 
GAEKGA (III, JJJ)  -  0.0 
180  CONTINUE 

190  CONTINUE 


....determine  the  EKGA  array . 

DO  220  I  -  1,6 
DO  215  J  -  1,6 
DO  210  K  -  1,6 

EKGA(I,J)  -  EKGA(I,J)  +  EKPR (I, K) *GAMMA (lEL, K, J) 
210  CONTINUE 

215  CONTINUE 

220  CONTINUE 

....determine  the  GAEKGA  array . 

DO  240  1-1,6 
DO  235  J  -  1,6 
DO  230  K  -  1,6 

GAEKGA(I,J)  -  GAEKGA(I,J) 

&  +GAMMA(IEL,K, I) *EKGA(K,  J) 

230  CONTINUE 

235  CONTINUE 

240  CONTINUE 
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C 


....copy  the  GAEKGA  array  into  the  EK  array . 

DO  260  1-1.6 
DO  250  J  -  1,6 

EK(IEL, I,J)  -  GAEKGA (I, J) 

250  CONTINUE 

260  CONTINUE 

270  CONTINUE 

....construct  the  GK  matrix . 

DO  300  lEL  -  1,  NEL 
II  -  3*{IEL-1) 

DO  290  J  -  1,  6 
JJ  -  II  +  J 

DO  280  K  -  1,  6 
KK  -  II  +  K 
GK(JJ,KK)  -  GK(JJ,KK) 

&  +EK(IEL, J,K)*{BASE(IEL)**NPOW) 

280  CONTINUE 

290  CONTINUE 

300  CONTINUE 

RETURN 

END 

************** ********************************************* ********** 


SUBROUTINE  FORCE  VECTOR  (NEL, NDOF, ELEN, ALPHA, BETA, FA, F) 


This  subroutine  is  used  to  construct  the  force  vector  for  the 
FEM  problem  specified. 


....declare  the  variables . 

INTEGER  NEL, NDOF, i, II, 12, 13, PI, P4 

PARAMETER (Pl-32,P4-99) 

REAL  ELEN, ALPHA (PI ), BETA (PI) , FA, F(P4) 

. . . .form  the  F-vector . 

F(l)  -  (ELEN/2.0)  *  (-COS (BETA (1) ) ) 

F(2)  -  (ELEN/2.0)  *  (COS (ALPHA (1) ) ) 

F(3)  -  (ELEN**2.0) /12.0 

DO  100  i-2,NEL 

11  -  (i-l)*3  +  1 

12  -  (i-l)*3  +  2 

13  -  (i-l)*3  +  3 

F(I1)  -  (ELEN/2.0)* (-COS (BETA (NEL))) 

&  + (ELEN/2.0)* (-COS (BETA (NEL-1))) 

F(I2)  -  (ELEN/2.0) * (COS (ALPHA (NEL) ) ) 

&  + (ELEN/2.0)* (COS (ALPHA(NEL-l))) 

F(I3)  -  0.0 
100  CONTINUE 

F(NDOF-2)  -  (ELEN/2.0)* (-COS (BETA (NEL) ) ) 
F(NDOF-l)  -  (ELEN/2.0) * (COS (ALPHA (NEL) ) ) 
F(NDOF)  -  -(ELEN**2.0) /12.0 

....scale  the  F-vector  by  FA . 

DO  200  i-l,NDOF 
F(i)  -  FA*F(i) 
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200  CONTINUE 


RETURN 

END 


*  * 
SUBROUTINE  BNDARY  (NDOF, GK, CLAN, FX, FY, FM, F, BXl , BYl , BMl , BX2, 

&  BY2,BM2) 

This  subroutine  is  used  to  in^ose  the  boundary  conditions  upon 
the  global  stiffness  matrix  and  force  vector. 

INTEGER  NDOF, BXl , BYl, BMl , BX2, 8X2, BM2, CLAN, i,N,Il,l2,l3,P4 
PARAMETER (P4-99) 

REAL  GK(P4,P4) ,FX,FY,FM,F(P4) 

....invoke  the  essential  boundary  conditions . 

IF  (BXl  .EQ.  1)  THEN 

Cy^LL  1MP0SE_BC  (NDOF, GK, 1,F) 

ENDIF 

IF  (BYl  .EQ.  1)  THEN 

CALL  IMPOSE_BC  (NDOF, GK, 2, F) 

ENDIF 


IF  (BMl  .EQ.  1)  THEN 

CALL  IMPOSE_BC  (NDOF, GK, 3, F) 

ENDIF 

IF  (BX2  .EQ.  1)  THEN 
N-NDOF-2 

CALL  IMPOSE_BC  (NDOF, GK, N, F) 

ENDIF 

IF  (BY2  .EQ.  1)  THEN 
N-NDOF-1 

CALL  IMPOSE_BC  (NDOF, GK, N, F) 

ENDIF 

IF  (BM2  .EQ.  1)  THEN 

CALL  IMPOSE_BC  (NDOF, GK, NDOF, F) 

ENDIF 

....add  the  concentrated  load  to  the  force  vector . 

11-  (CLAN-1) *3+1 

12-  (CLAN-1) *3+2 

13-  (CLAN-1) *3+3 

C  > 

F(11)-F(I1)+FX 
F(I2)-F(I2)+FY 

F(I3)-F(I3)+FM  • 

C 

RETURN 

END 


*  * 
SUBROUTINE  IMPOSE_BC  (NDOF, GK, N, F) 


This  subroutine  is  used  to  do  the  redundant  leg  work  of  inpos- 
ing  the  boundary  conditions. 
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C  ....declare  the  variables . 

INTEGER  ND0F,N,i,P4 
PARAMETER (P4-99) 

REAL  GK(P4,P4) ,F(P4) 

....impose  the  boundary  condition  on  the  GK  and  F  arrays . 

DO  100  i-l,NDOF 
GK(N,i)  -  0.0 
100  CONTINUE 

GK(N,N)  -  1.0 
F(N)  -  0.0 

RETURN 
END 

A******************************************************************** 

* 

SUBROUTINE  UPSCALE (NDOF, GK, F, BK, BF) 


This  subroutine  is  used  to  change  the  stiffness  matrix  &  force 
vector  from  single  precision  to  double  precision  in  order  to 
solve  the  linear  system  of  equations  in  double  precision. 


....declare  the  variables . 

INTEGER  NDOF,i, j,P4 
PARAMETER  (P4-99) 

REAL  GK(P4,P4) ,F(P4) 

REAL*8  BK(P4,P4) ,BF(P4) 

....generate  the  doubleprecision  compliments  of  GK  and  F . 

DO  110  i-l,NDOF 
DO  100  j-l,NDOF 

BK(i, j)  -  GK(i, j) 

100  CONTINUE 

BF(i)  -  F(i) 

110  CONTINUE 

RETURN 

END 

******************************************************************** 

* 

SUBROUTINE  DOWNSCALE (NDOF, BU, U) 


This  subroutine  is  used  to  do  down  scale  the  double  precision 
solution  of  the  linear  system  of  equations  back  to  single  pre¬ 
cision.  DOT  could  have  problems  with  double  precision  numbers! 


....declare  the  variables . 

INTEGER  NDOF,i,P4 
PARAMETER  (P4-99) 

REAL  U(P4) 

REAL* 8  BU(P4) 

....generate  the  doubleprecision  compliments  of  GK  and  F . 

DO  100  i-l,NDOF 
U(i)  -  BU(i) 

100  CONTINUE 
C 

RETURN 

END 

********************************************************************** 
*  * 
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SUBROUTINE  STRESS (X, Y, ALPHA, BETA, U, NEL, ELEN, YOUNG, H, SIGMA_T) 


This  subroutine  confutes  the  stress  at  each  nodal  point . 


. . . .declarations . 

INTEGER  NEL,NSNP,NDOF, SWITCH, i, II, 12, 13, 14, 15, 16, 17, 18, 19, 
&  P1,P2,P4 

PARAMETER (Pl-32, P2-33, P4-99) 

REAL  ELEN, H, X (P2 ) , Y {P2 ) , ALPHA (PI ) , BETA (PI ) , YOUNG, 

&  Kl,K2,CAl,CA2,CBl,CB2,v,vl,v2, 

&  elxdisp (P2) , elydisp (P2) ,  ELEN_f (P2 ) , DISPLEN (P2) , 

&  U(P4) ,SIGMA_N(P4) ,SIGMA_B(P4) , SIGMA_T(P4) 

....determine  the  constants . 

Kl-6.0/ (ELEN**2.0) 

K2-2.0/ (ELEN) 

NSNP  -  (NEL  +  1) 

NDOF  -  NSNP *3 

....determine  the  bending  stresses . 

DO  100  i-2,NEL 

11- (i-2)*3+l 

12- (i-2) *3+2 

13- (i-2)*3+3 

14- (i-l)*3+l 

15- (i-l)*3+2 

16- (i-l) *3+3 

17- (i-0) *3+1 

18- (i-0)*3+2 

19- (i-0)*3+3 
C 

CBl-  COS (BETA (i-1 ) ) 

CAl-  COS (ALPHA (i-1 ) ) 

CB2-  COS (BETA (i)) 

CA2-  COS (ALPHA (i ) ) 

C 

v2  -  Kl* (U(I4) -U(Il) ) *CB1 

&  +K1* (U(I2) -U(I5) ) *CA1 

&  +K2* (U(I3)+2.0*U(I6) ) 

C 

vl  -  Kl* (U(I4)-U(I7) )*CB2 

&  +K1*(U(I8)-U(I5))*CA2 

£  -K2* (U(I6)*2.0+U(I9) ) 

C 

IF  (ABS(vl)  .GE.  ABS(v2))  THEN 

V"Vl 

ELSE 

v-v2 

ENDIF 

C 

SIGMA_B(i)  -  YOUNG* (H/2.0) *v 
C 

100  CONTINUE 
C 

V  -  Kl*(U(l)-U(4))*COS(BETA(l)) 

£  +K1*(U(5)-U(2))*C0S(ALPHA(1)) 

£  -K2* (U(3) *2.0+U(6) ) 

C 

SIGMA_B(1)  -  YOUNG* (H/2.0 )*v 
C 

V  -  Kl* (U(NDOF-2)-U(NDOF-5) )*COS (BETA (NEL) ) 
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tr 


m 


r 


« 


c 


& 

& 


c 

c 

c 


c 


c 


c 


c 


300 

C 


C 


& 


c 


350 


+K1* (U (NDOF-4) -U (NDOF-1) ) *COS (ALPHA (NEL) ) 
+K2* (U(NDOF-3)+2.0*U(NDOF) ) 

SIGMA_B(NEL+1)  -  YOUNG* (H/2 . 0) *v 

....determine  the  normal  stresses . 

SWITCH  -  1 


IF  (SWITCH  .EQ.  1)  THEN 
DO  300  i-2,NEL 

11  -  (NEL-2)*3+l 

12  -  (NEL-2)*3+2 

13  -  (NEL-2)*3+3 

14  -  (NEL-1)*3+1 

15  -  (NEL-l)*3+2 

16  -  (NEL-l)*3+3 

17  -  (NEL-0)*3+1 

18  -  (NEL-0)*3+2 

19  -  (NEL-0)*3+3 

CAl-  COS (ALPHA (NEL-1 ) ) 

CBl-  COS (BETA (NEL-1 ) ) 

CA2-  COS (ALPHA (NEL) ) 

CB2-  COS (BETA (NEL) ) 

v2  -  (U(I4) -U(Il)  )  *CUV1  +  (U(I5) -U(I2)  )  *CB1 
vl  -  (U(I7)-U(14) )*CA2  +  (U(I8) -U(I5) ) *CB2 

IF  (ABS(vl)  .GE.  ABS(v2))  THEN 
v»vl 
ELSE 
v-v2 
ENDIF 

SIGMA_N(i)  -  (YOUNG/ELEN) *v 
CONTINUE 


V-  (U(4)-U(l))*COS(ALPHA(l))  +  (U (5) -U (2) ) *COS (BETA (1) ) 
SIGMA_N(1)  -  (YOUNG/ELEN) *v 

V  -  (U(NDOF-2)-U(NDOF-5) )*COS (ALPHA (NEL) )  + 

(U (NDOF-1 ) -U (NDOF-4 ) ) *COS (BETA (NEL) ) 

SIGMA  N(NEL+1)  -  (YOUNG/ELEN) *v 


ELSE 

DO  350  i-l,NEL+l 
SIGMA_N(i)  -  0.0 
CONTINUE 
ENDIF 


C 

C  ....determine  the  total  stresses  at  each  node.... 

DO  400  i-l,NEL+l 

SIGMA_T(i)  -  ABS(SIGMA_B(i))  +  ABS (SIGMA_N (i) ) 
400  CONTINUE 
C 


RETURN 

END 


* 

c 


SUBROUTINE  ARCH  OUTPUT 


* 
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This  subroutine  formats  the  final  results  and  output  of  the 
optimization  problem  and  stores  it  in  a  file  named  ARCH_OUT.DAT 


. . . .declare  variables. . . 

INCLUDE  'ARCH_COM.FOR' 

REAL  VOL, VOLUME 

....open  output  file  and  write  header . 

OPEN  (9,  FILE-'  ARCH_OUT .  DAT'  ,  STATUS-'  UNKNOVfN'  ) 

WRITE (9, 100)  LABEL 

WRITE (9, 100)  '  OPTIMIZATION  SOLUTION' 

WRITE(9,105)  '  - 

4 - • 

100  FORMAT (/5X, A) 

105  FORMAT {5X, A) 

. . . .section  "A" . 

WRITE (9, 100)  '  A)  Problem  Parameters:' 

WRITE (9, 110)  '  Arch  Angle  ANGLE,  '  Youngs  Modulus :', YOUNG 

WRITE{9,110)  '  Arch  Radius:',  RADIUS,  '  Yield  Strength: ', YIELD 

WRITE (9, 115)  '  Arch  Height:',  H,  '  No  of  Elements :', NEL 

110  FORMAT{8X, A,F12.3,T38,  A,F12.1) 

115  F0RMAT(8X,A,F12.3,T38,A,  110) 

....  section  ''B” . 

WRITE (9, 100)  '  B)  Derived  Constants:' 

WRITE (9, 120)  '  No  of  System  Nodal  Points ...', NSNP 
WRITE(9,120)  '  No  of  Degrees  of  Freedom ',NDOF 


WR1TE(9,125)  '  Length  per  Element . ',ELEN 

WRITE (9, 125)  '  Phi  Angle  per  Element . ',PHIANG 

WRITE (9, 120)  '  Number  of  Iterations . ', ITERATE 


120  FORMAT (8X, A, 16) 

125  FORMAT(8X,A,F12.4) 

....section  "C” . 

WRITE (9, 100)  '  C)  Structure  Loading:' 


WRITE  (9, 125)  'FX . '  ,  FX 

WRITE(9,125)  'FY . ',FY 

WRITE  (9, 125)  'FM . '  ,  FM 

WRITE  (9, 125)  'FA . ',FA 


....section  "D" . 

WRITE (9, 100)  '  D)  Elemental  Dimensions  and  Stress  Distribution:' 
WRITE(9,210)  'Element' , 'Height' , 'Base' , 'Length' , 'Volume' 

210  FORMAT(8X,A,T19,A,T34,A,T50,A,T62,A) 

220  FORMAT(8X,I4,T17,F10.5,T32,F10.5,T48,F8.5,T60,F8.5) 

VOLUME  -  0.0 

DO  300  i-l,NEL 

VOL  -  H*ELEN*BASE(i) 

WRITE(9,220)  i, H, BASE (i ), ELEN, VOL 
VOLUME  -  VOLUME  +  VOL 
300  CONTINUE 

. . . .section  "E" . 

WRITE (9, 100)  '  E)  Objective  Function:' 

WRITE(9,310)  '  Total  structure  Volume: ', VOLUME 
310  FORMAT(8X,A,F12.6) 
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WRITE(9,330)  'Node' , 'Stress' 

DO  320  i-l,NSNP 

WRITE (9,*)  i,SIGMA_T{i) 

320  CONTINUE 

330  FORMAT(8X,A,T19,A) 

....section  "F" . 

WRITE (9, 100)  '  F)  Boundary  Conditions;' 

WRITE (9, 410)  'Node' , 'X-Displ' , 'Y-Displ' ,  'Slope' 
WRITE (9, 430)  l,BXl,Byi,BMl 
WRITE(9,430)  NEL+1,BX2,BY2,BM2 

....  section  "G" . 

WRITE (9, 100)  '  G)  Solution  Vector:' 

WRITE (9, 410)  'Node' , 'X-Displ' , 'Y-Displ' , 'Slope' 
DO  400  i-l,NSNP 

11- (i-l) *3+1 

12- (i-l)*3+2 

13- (i-l)*3+3 

WRITE(9,420)  i,U(Il) ,U(I2) ,U(I3) 

400  CONTINUE 

410  FORMAT(T9,A,T17,A, T31,A,T46,A) 

420  FORMAT(7X,I5,3E14.6) 

430  FORMAT (7X, I5,T20, I4,T34, I4,T48, 14) 

RETURN 

END 


oo  o  o  ooooooooo 


ARCH  COMMON 


. . . .definitions . 

PI . The  maximum  number  of  elements 

P2 . The  maximum  number  of  global  nodal  points 

P3 . The  maximum  number  of  design  constraints 

P4 . The  maximum  number  of  degrees  of  freedom 

....declare  the  variables . 

INTEGER  NEL, NCON, NSNP, NDOF, METHOD, MINMAX,  INFO, IPRINT, IWK (400 
£  NRWK,NRIVnC, IPRM (20) , COUNT, OPTDCS, ITERATE, PRCSN, CLAN, 

£  BX1,BY1,BM1,BX2,BY2,BM2,P1,P2,P3,P4 


PARAMETER (P1-32,P2-33,P3-96,P4-99) 


REAL  ANGLE, RADIUS, ELEN, H,X (P2) ,Y{P2) , ALPHA{P1) ,BETA(P1) , 

£  YOUNG, YIELD, WK(27000) ,RPRM(20)  ,OBJ,G(P3) , 

£  DVIBG, DV1L0,DV1UP, BASE (PI ), BASEL (PI) , BASED (PI) , 

£  FA, FX, FY, FM, U (P4 ), SIGMA  T(P4) 


C 


£ 

£ 

£ 

£ 

£ 


. . . .make  in  common . 

COMMON  NEL, NCON, NSNP,NDOF, METHOD, MINMAX, INFO, IPRINT, IWK, 
NRWK,NRIWK, IPRM, COUNT, OPTDCS, ITERATE,  PRCSN, CLAN, 
BXl ,  BYl ,  BMl ,  BX2 ,  BY2 ,  BM2 , 

ANGLE, RADIUS, ELEN, H, X, Y, ALPHA, BETA, YOUNG, YIELD, 
WK, RPRM, OBJ, G, DVIBG, DVILO, DVIUP, BASE, BASEL, BASED, 
FA, FX,FY,FM,U, SIGMA  T 
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